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I.   INTRODUCTION 

A.   THE  MAIN  IDEA 

Image  processing  by  nonoptical  means  has  received  exten- 
sive attention  in  the  last  few  years.   Several  books  and  many 
papers  have  been  published  that  have  established  nonoptical 
image  processing  as  a  viable  area  of  research.   A  large 
portion  of  this  research  emphasizes  the  linear  processing 
of  images  for  two  main  reasons:   1)  Many  image  processing 
tasks  are  linear  in  nature.   These  tasks  include  image  enhance- 
ment, image  restoration,  picture  coding,  linear  pattern 
recognition,  and  TV  bandwidth  reduction.   2)  There  are  many 
known  linear  techniques  that  may  be  brought  to  bear  in  the 
treatment  of  linear  image  processing.   These  techniques  include 
transform  theory,  matrix  theory,  filtering,  signal  modeling, 
etc.   Several  common  operations  involved  in  image  processing 
include  transfer  function  concepts,  partial  difference  (recur- 
sive) equations,  and  convolution  summations.   For  example, 
Vander  Lugt  [Refs.  1,2]  has  presented  an  extensive  development 
of  linear  optics  based  on  transfer  functions.   The  transfer 
functions  relate  the  two-dimensional  Fourier  transform  of  an 
output  image  to  that  of  the  input  image.   Complex  optical 
systems  are  easily  described  by  combinations  of  transfer 
functions  that  correspond  to  individual  components  of  the 
optical  system. 


Partial  difference  equations  are  used  by  Habibi  [Ref.  3] 
to  describe  a  model  for  estimating  images  corrupted  by  noise. 
The  model  corresponds  to  a  two-dimensional  extension  of  Kalman 
filters.   Convolution  summations  are  discussed  by  Fryer  and 
Richmond  [Ref.  4]  in  work  that  involves  simplifying  a  two- 
dimensional  filter  to  a  single  dimensional  filter. 

The  time-discrete  state-space  model  offers  great 
utility  in  the  formulation  and  analysis  of  linear  sys- 
tems.  Linear  systems  that  are  described  by  transfer  functions, 
difference  equations,  or  convolution  summations  are  formulated 
into  a  state-space  representation.   Once  formulated,  many 
known  techniques  may  be  applied  to  systematically  analyze 
the  model.   Consequently,  the  state  space  model  is  a  general 
and  powerful  tool  that  is  used  to  unify  the  research  and  the 
study  of  time-discrete  linear  systems. 

This  thesis  develops  the  discrete  model  of  Roesser  [Ref. 
5]  for  linear  image  processing  which  closely  parallels 
the  well-known  state  space  model  for  time-discrete  systems. 
Because  it  is  parallel,  many  of  the  concepts  that  are  known 
for  the  temporal  model  may  be  carried  over  to  the  spatial 
model.   This  is  done  by  generalizing  from  a  single  coordinate 
in  time  to  two  coordinates  in  space.   The  spatial  model  will 
hopefully  have  some  of  the  same  utility  for  the  study  of  two- 
dimensional  linear  systems  as  the  temporal  model  for  one- 
dimensional  linear  systems  [Ref.  3].   However,  not  all  of  the 
properties  of  one-dimensional  systems  carry  over  into  the  multi- 
dimensional case. 
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One  of  the  fundamental  problems  involved  with  recursive  2- 
Dimensional  systems  is  that  the  order  of  the  system  (recursive 
memory)  is  not  the  same  as  the  number  of  initial  conditions 
(boundary  conditions) .   In  one-dimensional  systems  these  are  the 
same.   Temporal  systems  are  inherently  nonanticipatory  and  are 
often  treated  as  such  for  the  sake  of  physical  realizability 
in  real  time;  whereas  spatial  systems  do  not  have  causality  which 
is  an  inherent  limitation.   That  is,  an  image  processor  may  have 
right  to  left  dependency  as  well  as  left  to  right  dependency. 
Finally  it  is  noted  that  stability  criteria  in  one-dimensional 
recursive  systems  become  much  more  difficult  when  carried  over 
to  the  multidimensional  case. 

Causality  is  built  into  the  temporal  state-space  model  if 
an  initial  state  is  assumed  to  be  fully  specified.   In  order  to 
establish  a  close  parallel  for  the  spatial  model,  the  same 
built-in  causality  will  be  intentionally  assumed,  despite  the 
fact  that  causality  is  not  necessary  for  physical  realizability 
in  real  space.   Such  an  image  processor  is  said  to  be  unilateral. 
If  the  constraint  of  causality  is  removed,  then  the  image  proces- 
sor is  said  to  be  bilateral  [Ref.  5].   Concepts  that  are 
developed  in  this  thesis  for  the  latter  case  are: 

1)  Formulation  of  the  state  space  model  of  Roesser.  [Ref.  5] 

2)  The  definition  of  state  transition  matrix. 

3)  A  resulting  computer  program  based  on  the  above  model. 

4)  An  investigation  of  the  class  of  2-Dimensional  transfer 
functions  defined  by  this  model. 

5)  Derivation  of  a  general  response  formula. 


6)  Extension  of  Roesser's  model  of  state  variable  equations 
to  encompass  a  larger  class  of  transfer  functions. 

7)  Adaptation  of  the  1-D  "SSPACK"  program  to  produce  2-D  data 


B.   STATE  SPACE  REPRESENTATION 

Toward  the  end  of  the  19  50s,  the  concept  of  representing 
a  discrete  system  by  a  set  of  first-order  difference  equations 
became  a  standard  tool  of  the  research  engineer.   These  tech- 
niques have  since  become  generally  known  as  state-space  repre- 
sentations.  Such  representations  have  become  increasingly 
important  during  the  intervening  years  because^ they  often  allow 
one  to  carry  out  a  meaningful  system  design  entirely  in  the 
discrete-time  domain  (in  comparison  to  popular  Z-transform 
methods).   That  this  is  important  follows  basically  from  these 
factors : 

1.  The  system  may  be  nonlinear  so  that  transformation 
methods  are  not  directly  applicable. 

2.  Time-domain  concepts  often  give  one  a  better  insight 
into  the  analysis  and  synthesis  of  the  system  (fre- 
quently with  the  aid  of  a  digital  computer) . 

3.  Cases  in  which  the  initial  conditions  are  non-zero  may 
be  handled  straightforwardly. 

A  state  space  representation  of  a  system  differs  from  the 
conventional  representation.   In  a  conventional  representation 
only  the  relationships  between  the  input  and  output  signals  need 
be  known.   On  the  other  hand,  the  state-space  representation 
gives  a  total  description  of  both  the  internal  as  well  as  the 
external  signals  of  a  system. 
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C.   STATE- VARIABLE  REALIZATIONS — THE  CONCEPT  OF  STATE  — 
In  1-D  linear  systems  theory  and  control  theory,  the 
concept  of  a  filter  state  has  played  an  important  role. 
Basically  the  filter  state  at  any  point  in  time  contains 
all  the  information  necessary  to  compute  the  remainder  of  the 
filter  output  signal,  given  the  input  signal.   One  dimensional 
single-input,  single-output  filter  realizations  based  on  a 
state  variable  model  can  be  written  in  the  form: 

x(k+l)   =   Ax(k)  +  Bu(k)  (I. la) 

y(k)   =   Cx(k)  +  Du(k)  (I. lb) 

This  form  relates  the  input  u(k)  and  the  output  y(k)  through 
a  state  vector  x(k) .   The  state  vector  evolves  in  time  accord- 
ing to  equation  (I. la) .   The  matrices  A,  B,  and  C  and  1  xl  matrix 
D  govern  the  exact  form  of  the  input-output  relationship. 
(In  general  these  matrices  may  vary  with  the  index  (k)  and 
the  input  and  output  signals  may  be  vectors  as  well.)   Quite 
often  the  components  of  the  state  vector  are  taken  to  be  the 
constants  of  the  Z    delay  operators  in  a  flowgraph  represen- 
tation of  the  1-D  filter. 

A  classic  problem  in  state-variable  theory  representation  is  to 
find  the  matrices  A,  B,Cand  D  which  will  realize  a  particular  system 
function  H(z)  with  a  minimum  number  of  state  variables.   A 
similar  approach  may  be  taken  to  develop  a  2-D  state-variable 
model . 
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A  2-D  discrete  system  may  be  defined  as  a  mathematical 
abstraction  which  utilizes  three  types  of  variables  to  repre- 
sent or  model  the  dynamics  of  a  discrete-time  process.   The 
three  variables  are  called  the  input,  the  output,  and  the 
state  variable.   The  input  variables  u(i,j),  serve  as  external 
forces  which  influence  the  dynamics  or  motion  of  the  system. 
The  output  variables  y(i,j)  are  the  characteristic  variables 
which  are  directly  observable  (measurable)  by  the  external 
observer.   The  state  variables  x(i,j)  characterize  the  internal 
dynamics  of  the  system  and  are  to  be  selected  according  to  the 
following  rule. 

These  variables  are  formulated  in  such  a  manner  that,  if 
one  knows  the  values  of  the  present  state  variables  x(i,j)  along 
with  the  values  of  the  input  variables  u(i,j)  then  the  output 
variables  y(i,j>  and  the  next  state  variables  x(i,j)  are  com- 
pletely determined.   Moreover,  the  number  of  state  variables 
used  in  a  state-space  representation  must  be  minimized.   A 
state-space  representation  may  be  visualized  in  block  diagram 
form,  as  shown  below. 


u,ci> 

IL(U- 


.NTERNA  L 

STATES 

x,(U)   ,   W)...*^u 


Y.dJ ) 

Y,C  l.J  ) 
Y,(U) 


Figure    1.1 
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In  Figure  1. 1 ,  m- inputs ,  p-outputs and  n-state  variables  are  repre- 
sented.  However,  we  will  be  mainly  interested  in  those  systems  which 
have  one  input  (m  =  1)  and  one  output  (p  =  1) .   It  is  important 
to  note  that  the  input  and  output  variables  appear  external  to 
the  system,  while  the  state  variables  are  generally  internal. 

The  different  input  variables  will  be  represented  by  the 
input  vector  u(i,j)  where, 


u(i, j) 


u1(i,  j) 
u2 (i,  j) 


u  ( i  ,  j  ) 

m 


the  output  vector  y(i,j)  where/ 


y(i. j) 


Y1(i.j) 

y2Cif  j) 


yp(i,j) 


and  the  state  vector  x(i,j)  where, 


x(i,  j) 


x2  (i,  j) 


xn(i,j) 


For  a  given  process  the  state  space  representation  is  not  unique 
However  all  such  representations  have  one  characteristic  in 
common  for  a  given  system,  namely  the  number  of  elements  n  is 
referred  to  as  the  order  of  the  system. 
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II.   ROESSER'S  STATE-SPACE  MODEL 

A.   THE  FRAMEWORK 

An  image  is  a  generalization  of  a  temporal  signal,  in  that 
it  is  defined  over  two  spatial  dimensions  instead  of  a  single 
temporal  dimension.   Consequently,  two  space  coordinates  i  and 
j  take  the  place  of  time,  t.   Also,  two-state  sets  are  intro- 
duced to  replace  the  single-state  set.   The  following  defini- 
tions are  made  by  the  model: 

i     An  integer-valued  vertical  coordinate; 

j     An  integer-valued  horizontal  coordinate; 

{R}   A  set  of  n^  real  vectors  which  convey  information 
horizontally ; 

{S}   A  set  of  n~  real  vectors  which  convey  information 
vertically; 

{u}   A  set  of  m  real  vectors  that  act  as  inputs; 

{y}   A  set  of  p  real  vectors  that  act  as  outputs. 

A  specific  image  processor  is  then  defined  as  6-tuple 

<{R},{S},{u}, {y},f ,g>  , 

where  f  is  the  next  state  function: 

f:   {{R},(S},{u}   -   {{R},(S}} 
and  y  is  the  output  function 

g:   {{R},{S},{u}}   -»   {y}  . 
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Now  since  f  and  g  are  to  be  linear  functions,  they  may  be 
represented  by  the  following  matrix  equations: 

R(i+l,j)   =   AJ_R(i,j)  +A2S(i,j)  +  B1u(i,j) 

S(i,j+1)   =   A3R(i,j)  +  A4S(i,j)  +  B2u(i,j)      (II. 1) 

y(i,j)     =   C1R(i,j)  +  C2S(i,j)  +  Du(i,j)   i,j  >  0 

Ai ,  A„ ,  A.,,  A.,  B,,  B„ ,  C, ,  C2 ,  D  are  matrices  of  appropriate 
dimensions.   Boundary  conditions  R(0,j)  and  S(i,0)  and  also  the 
input  u(i,j)  are  externally  specified.   In  the  next  section  a 
computational  rule  is  obtained  that  uniquely  determines  the 
states  R(i,j)  and  S(i,j)  and  also  the  output  y(i,j)  (for  i,j  >  0 
from  the  boundary  conditions  (such  as  all  zero).   The  equations 
produce  a  set  of  output  vectors  from  the  input  vectors. 

This  formulation  is  general  so  that  any  discrete  linear 
image  process  may  be  so  represented.  Notation  is  condensed 
somewhat  by  introducing  the  following  matrices  and  vectors: 


A   = 


T(i,j) 


Al   A2 
A3   A4 


B   = 


R(i,j) 
S(i,j) 


B. 


B. 


T' (i,j) 


c  =  [c1  c2 


R(i+l,j 
s(i,  j  +  1. 


T' (i,j)   =   AT(i,j)  +  Bu(i, j) 

y(i,j)   =   CT(i,j)  +  Du(i,j) 

B.   GENERAL  RESPONSE  FORMULA 

A  state-transition  matrix  A  is  defined  as  follows 
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A   = 


Al   A2 

A3   A4 


Then  exponentiation  A     is  defined  as, 


A1^   =   A^V"1^  +  A^V'^1      (i,j)  >  (0,0) 

A0,0   =   I  ;  A~1,j   =   A1/_j   =   0    for   j  >  1,  i  >  1 


Examination  of  this  definition  bears  out  that  it  is  an  effective 
recursive  definition  of  A    for  integer  values  of  i  and  j  such  that 
either  i  >  0  or  j  >  0  or  (i,j)  =  (0,0).   It  parallels  the  defini- 
tion of  the  time-discrete  state-transition  matrix  A   =  A  A 

It  now  remains  to  be  shown  that  this  state  transition  matrix 
A  /J  may  be  used  in  expressions  for  the  response  of  the  model 
in  terms  of  the  inputs  and  boundary  conditions.   The  term 
boundary  conditions  is  used  here  to  refer  to  the  states  along 
the  edges  of  the  model.   Specifically,  the  set  of  boundary  condi- 
tions consist  of  R(0,j)  for  j  _>  0  and  S(i,0)  for  i  >_   0. 

C.   CHARACTERISTIC  FUNCTION  OF  A  MATRIX 

If  the  primary  inputs  and  outputs  are  dropped  in  the  model 
equations  (II. 1),  a  representation  arises  for  the  state  behavior 
of  the  system  having  the  form 


R(i+l,j)   =   A,R(i,j)  +  A0S(i,j 


S(i,j+1)   =   A3R(i,j)  +A4S(i,j) 


(II. 2) 


These  equations  are  useful  in  the  development  of  a  form  for  a 
two-dimensional  characteristic  matrix  of  A.   Operators  are 
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first  introduced  that  advance  a  particular  coordinate  of 
their  operand. 

Definition:   Let  E  be  an  operator  that  has  the  effect  of  ad- 
vancing the  vertical  coordinate  or  the  first  subscipt  of  the 
function  upon  which  it  is  operating.   Likewise,  let  F  be  an  opera- 
tor that  has  the  effect  of  advancing  the  horizontal  coordinate  or 
second  subscript  of  the  function  upon  which  it  is  operating. 
The  effect  of  these  operators  on  the  state  vectors  is: 

R(i+l,j)   =   ER(i,j) 
S(i,j+1)   =   FS(i,j) 

The  state  equations  can  be  rewritten  using  these  advance 
operators . 


(EI-Ax)R(i, j)  -  A2S(i,j)   =   0 


-A3R(i,j)  +  (FI-A4)S(i,j)   =   0 


These  equations  are  equivalently  represented  in  the  overall 
matrix  form. 


(EI-A-,)      -A, 


-A. 


(FI-A 


T(i,  j 


=   0 


The  above  equation  represents  a  system  of  homogeneous  equations  in 
the  elements  of  T(i,j).   If  the  system  is  to  have  a  non-trivial 
solution  for  T(i,j)  then  the  transformation  represented  by  the 
matrix  must  be  singular.   The  above  matrix  is  said  to  be  the  two- 
dimensonal  characteristic  matrix  of  the  partitioned  matrix  A,  where 
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A   = 


Al   A2 
A3   A4 


The  characteristic  matrix  of  A  is  denoted  cm (A)  and  may  be 
represented  as 

cm (A)   =   EI1,0  +  FI0/1  -  A 

where, 


.1,0 


and 


.0,1 


0  I  0 

o" ;  x 


Now  since  cm (A)  must  be  singular,  its  determinant  must  be  equal 
to  zero.   | cm (A) |  =  0.   If  E  and  F  are  replaced  in  the  above  by 
general  indeterminates  x  and  y  respectively,  the  result  is  an 
expression  called  the  two-dimensional  characteristic  equation 
for  A.   The  determinant  of  cm (A),  and  x  and  y  replacing  E  and 
F,  is  called  the  two-dimensional  characteristic  function  of  the 
matrix  and  is  denoted  by 


cm(A) 


f(x,y)   =   0 


f (x,y)  will  be  a  monic  polynomial  in  x  and  y  with  degree  n, 
in  x,  and  degree  n~  in  y,  where  n,  is  the  dimension  of  R  and 
n„  is  the  dimension  of  S.   f(x,y)  has  the  form 


f (x,y) 


(0,0)  <  (i,  j)  <  (n-,,n2) 


a .  .x  yJ 


where  a.  •  denotes  elements  of  A  and  a       =  1. 
lfD  nl'n2 
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D.   CIRCUIT  ELEMENTS  AND  THEIR  REALIZATION 

Let  us  consider  the  single  2-D  IIR  filter  transfer  function 
given  by: 


H(Z1,Z2) 


b00+b10Zl1+b01z21+bll2Ilz21+b21ZI2z21 
~     1=1     rI     -1  -1     -2  -1 — 
1~a102l  "a01Z2  "allZl  z2  "a21Zl  z2 


B(z,,z2) 

1-A(zlf z2) 


A  simple  block  diagram  for  H(z,,z  )  follows. 


uLCiJ) 


B  (2.,  20 


y(u; 


W(i,JJ 


A  U,,lO 


Figure  2.1 

The  input  signal  u(i.,j)  flows  through  a  filter  corresponding 
to  the  numerator  transfer  function  B(z,,z  ).   The  resulting 
signal  is  added  to  the  signal -w (i , j )  to  produce  the  output 
signal  y(i,j).   The  denominator  transfer  function  1-A(z,,z  ) 
is  realized  by  the  feedback  loop  containing  A(z,,z  ). 

Since  we  are  dealing  with  two  dimensions,  there  are  two 
fundamental  shift  operators  which  may  occur  along  a  signal  flow 
path,  the  horizontal  shift  operator  indicated  by  z,   and  the 
vertical  shift  indicated  by  z2   [we  shall  omit  from  consideration 
the  inverse  shift  operators  z-,  and  z2]  .   In  most  cases  of  prac- 
tical interest  they  can  be  eliminated  by  multiplying  both  the 
numerator  and  denominator  polynomials  of  H(z,,z2)  by  the  appro- 
priate powers  of  z,l  and  z?  .   Let  us  look  at  a  signal  flowgraph 
representing  the  numerator  polynomial: 
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B(Z1'Z2}   =   b00  +  b10Zl1  +  b01Z21  +  bllZl3"Z21  +  b21Zl2z21 


which  is  shown  in  Figure  2.2  below. 


Figure  2.2 


Note  the  chain  of  two  z,   operators  descending  on  the  left  and 
the  single  z2   operator  ascending  on  the  right.   The  nodes 
along  these  two  vertical  paths  are  connected  by  branches  with 


the  appropriate  gains.   If  we  label  the  nodes  in  both  z. 


-1 


chains  and  the  z~   chain  0,1,2  and  so  on,  from  the  top  down, 
the  ith  node  in  the  z~   chain  is  connected  to  the  jth  node 

in  the  z^1  chain  by  a  branch  with  a  gain  factor  of  b. .. 
2  J  3  lj 

Similarly  the  signal  flowgraph  for  the  polynomial  A(z,,z 
is  shown  in  Figure  2.3. 


\N(»0"T 


OWTfWl 


Figure    2 . 3 
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Since  there  is  no  aQQ  term,  there  is  no  direct  connection  be- 
tween the  input  and  output  nodes  of  this  signal  flowgraph.   Thus 
any  path  from  the  input  node  to  the  output  node  will  encounter 
at  least  one  z-,   or  z„   shift  operator. 

At  this  point  it  is  appropriate  to  discuss  realizations  for 
the  two  shift  operators  z,   and  z2  .   At  their  simplest  level, 
the  shift  operators  merely  select  the  "previous"  S-tuple  value 
in  the  horizontal  or  vertical  direction.   When  the  input  to  a 
z,   operator  is  the  S-tuple  u(i,j)  the  output  will  be  R(i-l,j). 
Similarly  for  a  z_   operator  the  output  will  be  S(i,j-1)  when 
the  input  is  R(i,j)  or  S(i,j).   Consequently  a  realization  of 
either  shift  operator  must  embody  the  appropriate  amount  of  memory 
to  retain  the  "previous"  S-tuple  in  the  appropriate  direction. 

Interestingly  enough,  in  the  more  general  case  where  the  numer- 
ator and  denominator  polynomials  are  considered  jointly,  the  state 
variable  realizations  based  on  conventional  signal  flowgraphs  may 
not  be  minimal  in  the  sense  that  the  transfer  furnction  can  be 
realized  with  fewer  coefficients.   Consider, 


H(z1,z2) 


b10Zl  +b01Z2  +bllZl  Z2 

: =i =i -l  -l 

1_a10Zl  "a01Z2   allZl  Z2 


(II. 3) 


The  corresponding  signal  flow  representation  is  shown  in  Figure 
2.4  below: 


yg^ 
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E.   ANALYSIS  OF  ROESSER'S  MODEL 

Recalling  from  page  14  the  equations  of  the  model  are: 

R(i+l,j)   =   AxR(i,j)  +  A2S(i,j)  +  B1u(i,j) 

S(i,j+1)   =   A3R(i,j)  +  A4S(i,j)  +  B2u(i,j) 

y(i,j)   =   C1R(i,j)  +  C2S(i,j)  +  Du(i,j) 

A,,  A~  ,  A_ ,  K.,    B,  ,  B2 /  C, ,  C~ ,  D  are  scalars  or  matrices  of 
appropriate  dimensions. 


R(i+1, j) 
S(i, j+1) 


Al   A2 
A3   A4 


R(i,j) 
S(i,j) 


+ 


B 


B. 


u(i,j)      (II. 4) 


y(i,j)   =   [c1  c2] 


R(i,j) 
S(i,j) 


+  Du(i,j) 


(II. 5) 


R(i+l,j)   =   A]_R(i,j)  +  A2S(i,j)  +  Bxu(i,j) 


S(i+l,j)   =   A-.R(i,j)  +  A.S(i,j)  +  B9u(i,j) 


And  taking  Z  transforms: 


z1R(z1,z2)   =   A1R(z1,z2)  +  Z2S(z1,z2)  +  B1u(z1,z2 


z2S(z1,z2)   =   A3R(z1,z2)  +  A4S(z1,z  )  +  B^fz^z  ) 


y(z1,z2)   =   [C±      C21 


R(z1# z2) 
S(z,,z  ) 


+  Du(zlfz  ) 


(II. 6 
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or 


Z1R(Z1,Z2)     "   AiR(z!'z2)    "   A2S(z1,z20       =      B1u(z1,z2) 
z2S(z1,z2)     -    A3R(z1,z2)     -    A4S(z1,z2)       =      B2u(z     ,z2) 


or 


R(z1,z2)  [z1  -A1J     -    A2S(z1,z2)       =      B1u(z1/z 


R(z,,z    )[-A7]     -     [z„-A]S (z, ,z_)    =      B0u(zn,z0] 


'l'"2 


l'"2 


'2~v"l'    2 


or 


Zi-Ai 


-A, 


-A, 


Z2-A 


R(zx,z2) 
S (z^z    ) 


B. 


B. 


u(z1,z2 


or 


zii  ° 

VJ~z2~ 


Al  i  A2 
V3  "•  Y4 


R(zlfz    ) 

sTz1,z2) 


B- 


3. 


u(zl'z2 


where    Z-    =    z, I    and    Z_    =    z_I,    and 


R(z1,z    ) 
S  (z1,z    ) 


z    !   0 

0      '"? 


A- 


B„ 


u (z1# z2) 


and    after   substitution    in   Equation    (II. 6) 


y(zl'z2}       =       CC 


1      C2] 


0 


*1  ' 

-  "  1 
0     ,  z 


I  "2 


1 

-1 

1 

Ai ;  A2 

Bi 

A3   J  A4 

_BJ 

u  (z^z- 


+   Du(z1/Z2) 
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or 


H(z1,z2) 


y (z1/ Z2) 
u(z1,z2) 


=   [C 


+  D 


c2] 


Zl  -  ° 
0  I  z0 


Al' 

A2 

-1 

Bl 

A3 

A4_ 

B2 

1 



1         



(II. 7) 


Th 


e  submatrix  Z,  is  simply  z,  times  an  identity  matrix  of  the 


appropriate  size.   Similarly  Z„  is  z„  times  an. identity  matrix. 
The  objective  of  the  state  variable  realization  procedure  is 
to  find  the  matrices  A,  B,  C,  and  D  which  yields  an  F(z ,,z  ) 
that  equals  or  approximates  a  desired  system  function  H(z,,z  ). 
In  essence,  the  equations  of  Roesser  represent  an  implementation 
for  which  a  design  algorithm  must  be  found.   One  choice  for 
the  state  variables  is  the  output  signals  from  the  shift 
operators. 

Thus  R(i,j)  is  a  vector  containing  the  output  signals 
from  the  z..   operators  and  S(i,j)  contains  the  output  signals 
from  the  z„   operators.   (Note  that  the  output  signal  of  a 
shift  operator  signal  path  is  not  necessarily  the  same  as  the 
nodal  signal  at  the  node  to  which  the  signal  path  points.) 
If  a  state  variable  corresponds  to  the  output  of  a  shift  operator, 
the  next  value  of  that  state  variable  must  correspond  to  the 
input  of  the  shift  operator.   To  obtain  the  submatrices  A,, 
A2 1    A-.,  A.  in  equations  of  Roesser,  we  write  the  input  signal 
of  each  shift  operator  in  terms  of  the  outputs  of  all  the 
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shift  operators,  taking  care  to  include  all  shift-free  paths 
from  ouput  to  input  (see  the  following  flowgraph) . 


U(U) 


YCU) 


Figure  2.5 


Expanding  the  form  of  Equation  II-7,  page  24,  yields 


H<VZ2 


[C1      C2] 


\ 

0 

Al 

1*2 

-1 

Bi 

—  — 

—  — 

—   — 

1  _     _ 

m~  «• 

0 

_ 

ZJ 

b 

1    A 
i     1 

i_B2_ 

+    D 


-1 


det   A 


adj    A 
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A 


-1 


[C1   C2] 


(z1-A1) (z2-A4)-A2A3 


JCi  J  _*2  . 


A. 


VAi 


3. 


3. 


[C1   C2] 


(Z2-A4)BX 


A2B2 


(Z2-A1) (A2-Z4)-A2A3 


A3B1 


Z1-A1)(Z2-A4)-A2A3 


(Z1-A1)B2 


(Z1-A1)(Z2-A4)-A2A3 


(Z1-A1)(Z2-A4)-A2A3 


[C1   C2] 


(Z2-A4)BX  +  A2B2 


Z2~A1) (Z2-A4) ~A2A3 


A3EX  +  (Z1-A1)B2 


(Z1-A1) (Z2-A4) ~A2A3 


or 


H  (Z1,Z2) 


C1(Z2~A4)B1  +  C1A2B2  +  C2A3B1  +  C2(Z1~A1)B2 


or 


H(z1/Z2) 


C1B1Z2-C1B1A4+C1A2B2+C2A3B1+C2B2Z1-C2B2A1 
Z1Z2-A4Z1    -    AlZ2-A2A3+AlA4 


(C1A2B2+C2A3B1-C2B2A1-C1B1A4)+(C2B2Z1+C1B1Z2) 


(AlA4-A2A3)     -    A4ZX   "A1Z2+Z1Z2 


Equating    equation    (II. 8)    with    (II. 3)    on   page    21    yields 


(II. 8) 
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(C1A2B2+C2A2B1-C2B2A1-C1B1A4)+C2B2Z1+C1B1Z2 

uyvw  -  a4zx  -  axz2  T  ^z2 


K  "I.!.  "I      i-  "I      "I 

b10Zl  +b01Z2  +bllZl  Z2 


i      "1      "I      "I  "I 
1  a10Zl  "a01Z2  "allZl  Z2 


For  this  example,  Z,  =  z.,  ,  Z„  =  z2,  all  of  the  coefficients  on 
the  left  hand  side  are  scalars.   Equation  terms  of  equal  powers 
of  z,  and  z~/ 


ClA2B2+C2A3BrC2B2ArClBlA4   =   bll   =   °  "  {11'9) 

C2B^   =   b1Q  (11.10) 

C1B1  =  boi  .(ii-ii) 

A1A4  "  A2A3   =  1  (11.12) 

(11.13) 
(11.14) 
(11.15) 


From  these  equations,  assuming  that  B,  =  B2  =  1,  it  follows  that 

Cl   =   b10 
C2   =   b0l 

Al  =  aoi 

A4   =   al0 
27 


A4 

aio 

Al 

= 

a0l 

an 

= 

-1 

From  Equation  (11-12)  : 


A1A4  "  A2A3   =  1      =      "all 


■A2A3   =   "all  "  A1A4 


A2A3   =    au  +  a10aQ1 


Let  A~  and  A_  take  on  particular  values  p  and  q  respectively, 

a3  =  q     A2  =  p 


or 


Pq   =   au  +  a10aQl  (11.16) 

From  Equation  (II-6)  : 


C1A2B2  +  C2A3B1  "  C2B2A1  "  C1B1A4   =   bll 


or, 


boiP  +  bioq  -  bioaoi  "  boiaio  "  bll  =  °      (I1-17 


Substituting  Equation  (11.16)  into  Equation  (11.17) 


an+aioaoi 

b01  q +  b10q  "  b10a0l  "  b0ia!0  "  bll   =   ° 


or 
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bioq2-(bioaoi+boiaio+bn)]   +  (boian+boiaioaoi)     =     °   •      (II-18) 


The  results  are  just  the  same  as  in  [Ref .  8] .   After  the  com- 
parison between  Roesser's  model  and  the  2-D  IIR  filter, 
described  by  Equation  (II. 3),  we  have: 


Al  "  a0l 

A2  =  p  ;   pq  =  axl  =  a^a^ 

A3  =  q;   b10q2-(b10a01+b0iai0+bll)q+(b0lfll+b0iai0a01 

A4  =  al0 

Cl  =  b01  (I1-19 

C2  =  b10 

Bl  =  1 

B2  =  1 

D  =  0 


=   0 


The  foregoing  equations  relate  the  coefficients  of  the  2-D 
transfer  function  to  the  terms  of  the  system  matrices  of  the 
Roesser  model,  Equation  (II. 1). 

Kung  et  al.  [Ref.  2]  have  shown  that  the  following  state 
variable  equations,  which  use  only  two  shift  operators,  will 
also  realize  H(z,,z  ).   For  the  foregoing  example, 
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R(i+l,j) 
S(i,j+1) 


Loi 


,  aio 


R(i,j) 
S(i,j) 


+ 


u(i,j) 


Y(i,j)      = 


[b10      b0l] 


R(i/j) 
S(i,j) 


or 


H(zrz2 


[bl0      b01] 


zi  aio 


-q 


-p 


Z2'a01 


-  -1 

-  -. 

1 

1 

-  — 

We  can  construct  a  signal  flowgraph  with  only  two  shift 
operators.   It  is  an  equivalent  figure  to  that  on  page  25. 

Kung  et  al.  [Ref.  2]  have  also  shown  that  state-variable 
realizations  of  the  form  of  the  equations  above  may  be 
generalized  for  any  system  function  H(z,,z„)  which  satisfies 
the  following  three  conditions: 


U(U) 


Y(U) 
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1)  The  constant  term  in  the  numerator,  fc>n0  =  0,  must  be 
zero. 

2)  The  largest  powers  of  z±    ,  in  the  numerator  and 
denominator  polynomials,  must  be  equal,  and 

3)  The  largest  powers  of  Z2   in  the  numerator  and 
denominator  polynomials  must  be  equal. 

There  is  one  potential  difficulty  with  state  variable  reali- 
zations of  this  type.   The  nonlinear  equations  defining  p  and 
q  may  result  in  complex  values  for  these  constants.   For 
example,  when  b, n  =  b_,  =  1,  b, ,  =  0,  a, n  =  a„,  =  2  and 
a1 .  =  1,  we  get  p  =  q*  =  2  ± j. 
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III.   THE  PROGRAM  OF  ROESSER'S  EQUATIONS 

WITH  SCALAR  COEFFICIENTS  (FIRST  ORDER) 

A.   AN  EXAMPLE 

For  a  4  x  4  data  field  the  S  and  R  matrices  are  indexed 
as  follows: 


j 

j 

1,1 

1,2         1,3 

1,4 

1,1 

1,2         1,3 

1,4 

2,1 

2,2         2,3 

2,4 

i 

2,1 

2,2         2,3 

2,4 

3,1 

3,2         3,3 

3,4 

1 

i 

3,1 

3,2         3,3 

3,4 

4,1 

4,2         4,3 

4,4 

4,1 

4,2         4,3 

4,4 

S   matrix 

R  matrix 

For  4x4  Matrices 
The  Initial  Conditions  are  given  by  the  values 

R(l,l) ,  R(2,l) ,  R(3,l) ,  R(4,l) 
S(l,l) ,  S(l,2) ,  S(l,3) ,  S(l,4) 

The  2-D  state  variable  equations  can  be  written  as: 


R(i+1, j) 
S(i,  j  +  1) 

yd,j) 


A-,R(i,j)  +  A0S(i,j)  +  B,u(i,j) 


=   A,R(i,j)  +  A.S(i,j)  +  B„u(i,j) 


[C1   C2] 


R(i,j) 

S(i, j) 


The  input  2-D  field  is  taken  to  be, 


u(i, j)   =   1,   for   i  =  j  =  1 
=   0,   otherwise. 
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The  output  data  field  is  indexed  as: 


2 

1,1  1,2  1,3  1,4 

2,1  2,2  2,3  2,4 

3,1  3,2  3,3  3,3 

4,1  4,2  4,3  4,4 


Y  output  matrix 

B.   THE  2-D  FOURIER  TRANSFORM 

The  2-D  discrete  Fourier  transform  Y(m,n)  gf  the  output 
y(i,j)  can  be  written  as, 


M-l  N-l         -J27T—   -j2tt— 
Y(m,n)   =    I         J      y(£,k)e     M  e     N 
£=0  k=0 


or    for   convenience, 


Y(m,n) 


M         N 

I         I    y(£,k)e 
£=1   k=l 


.„     U-lMm-l)          .       (k-1)  (n-l 
H2tt n -d2tt 


M 


N 


Y(m,n):       2-D   D.F.T.     {y(i,j)} 


M  xN:   The  dimension  of  the  given  data  y(£,k)  and  D.F.T 
Y(m,n)  also. 


y(£,k):   Given  data  (The  output  as  described  above). 

To  develop  the  D.F.T.  for  two-dimensional  signals  we 
consider  a  finite  area  sequence  y(£,k)  which  is  zero  outside 
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the  interval  0  <_  I    <_  M-l,  0  <_  k  <_   N-l,  i.e.,  it  is  of  area 
(M,N)  and  construct  the  periodic  sequence: 


yU,k)   =  y[((A))M(Ck))  )J 


The  original  sequence  y(£,k)  is  recovered  by  extracting  one 
period  of  y(£,k),  i.e., 


y(£,k)   =  yU,k)RM^NU,k) 


1,   0  <_   Z    <_   M-l,  0  <_.  k  <_   N-l 

M,  1m 

0,   otherwise  . 


We  then  define  the  discrete  Fourier  transform  of  y(£,k)  to 
correspond  to  the  Fourier  series  coefficients  of  y(£,k). 
However,  just  as  we  did  with  one-dimensional  sequences,  we 
will  maintain  the  duality  between  the  time  and  frequency 
domains  by  interpreting  the  D.F.T.  coefficients  to  also  be  a 
finite  2-D  sequence.   Thus  with  Y(m,n)  denoting  the  D.F.T. 
of  y(£,k) ,  we  can  write 


M-l    N-l  "J2^      ~^2TrTT 

Y(m,n)       =         I         I      y(£,k)e  e  R._  M(m,n: 

£^0    k=0  M'N 


or 


.     M-l    N-l  j2tt^      j2tt]^ 

y(£,k)       =     ±      I        I      Y(m,n)    e  M    e  H    RM      (£,k) 

m=0    y=0  ' 
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or, 


M    N  _j2TTU-l)(m-l)   _j2:;(k- 


Y(m,n)       I         I      y.  .(£,k) 
£=1  k=l    ': 


N         J        N 


As  an  example,  consider  the  case  for  M  =  N  =  5 
Given  2-D  Data  Sequence 


1,1  1,2  1,3  1,4  1,5 

2,1  2,2  2,3  2,4  2,5                       Matrix                  5x5 

yU,k)=              3,1  3,2  3,3  3,4  3,5                             M=5                     N=5 

4,1  4,2  4,3  4,4  4,5                 £=1,2,3,4,5      k=l,2,3,4,5 

5,1  5,2  5,3  5,4  5,5 


Then, 


y(l,l)  +  y(l, 2)  +  y(l, 3)    +  yd, 4)  +  yd, 5) 

+y(2,l)  +  y(2,2)  +  y(2,3)    +  y(2,4)  +  y(2,5) 

Y(l,l)    =         +y(3,l)  +  y(3,2)  +  y(3,3)   +  y(3,4)  +  y(3,5) 

10=1/11=1  +y(4,l)  +  y(4,2)  +  y(4,3)   +  y(4,4)  +  y(4,5) 

+y(5,l)  +  y(5,2)  +  y(5,3)    +  y(5,4)  +  y(5,5) 
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In  Appendix  A  we  give  a  listing  of  the  programs  that  have 
been  written  to  generate  y(i,j)  and  Y(m,n). 

C.   NUMERICAL  EXAMPLES 

Three  numerical  examples  which  depend  on  Equation  (11.19) 
are  used  to  demonstrate  the  program  in  Appendix  A. 
First  example; 

^(z"1  +Z"1) 
H(z.,z0)   = 


1  -  .  2Z-,   -  .  3z„ 


yields 


all   =   ° 


Al   =   a01   =   °'3 


A4   =   a1Q   =   0.2 


Cl   =   b10   =   °'5 


C2   =   bQl   =   0.5 


Bl   =  1 


B2   =   1 


D    =   0 


After  substitution  of  these  values  in  Eqs .  (11.16)  and  (II. 1! 
we  identify 


au   =   0    bQ0   =   0     bxl   =   0 


B- 


B. 


a0l 

— 

0.3 

p    = 

0 

.2 

q    = 

0 

.3 

aio 

= 

0.2 

boi 

= 

0.5 

b10 

= 

0.5 

1 

1 

0 

(III. la) 


Second  example: 

Proceeding  in  a  similar  way  with 


H(z1,z2) 


0.25z  1+0.3z21+0.2z11z2 
l-0.125z~1-0.2z21-0.1z~1z~1 


yields 


a .,   =   0.1 


b00   =   ° 


aQl   =   0.2 


bxl   =   0.2 


A, 


p   =   0.125 


=   0.83 


q  =   1 


or    =   0.15 


(III. lb 


a1Q   =   0.125 


C 


bQl   -   0.3 


b1Q   =   0.25 
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Bl  =  1 


B2   =   1 


D    =   0 


Third  example: 

Proceeding  in  a  similar  way  with 


0.2  5z  X+0.15z  1+0.72z11z~1 
K(z,  ,z9) 


1 '  2'  -l       -1       -1-1 

l-0.13  5z1  -0.2  5z2  -0.15z1  z2 


yields 


bQ0   =   0     bu   -   0.25 


all   =   °*15 


Al  =  a0l   =   °'25 

A2  =  p   =   0.1312 

A3  =  q   =   1.4  (III.lc) 

A4  =  al0   =   0.135 

ci  =  boi  =  °-15 

C2  =  blQ   =   0.25 

Bl  =  1 

B2  =  1 

D  =  0 


Zero  initial  conditions  were  assumed  for  all  examples.   The 
simulation  results  are  presented  in  Figures  3-1,  3-2  and  3-3 
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In  order  to  verify  the  correctness  of  the  output  produced 
by  Roesser,  the  2-D  D.F.T.   Y(m,n)   plots  for  these  examples 
were  compared  with  the  corresponding  |h(z,,z  ) |.   2-D  transfer 
function  plots  |H(z,,z  ) |  for  Examples  1,  2  and  3  are  shown 
in  Figs.  3-4a,b,  3-5a,b  and  3-6a,b  respectively.   The  listing 
of  a  program  used  to  generate  these  plots  can  be  found  in 
Appendix  B. 
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IV.   EXTENSION  OF  ROESSER'S  MODEL  TO  SECOND  AND  HIGHER  ORDERS 

A.   MINIMIZING  THE  NUMBER  OF  SHIFT  OPERATORS 

In  order  to  minimize  the  number  of  shift  operators  we  follow 
the  procedure  given  in  Kung  [Ref.  8].   Let  us  consider  the 
simple  2-D  IIR  filter  transfer  function  given  by 


H(z1,z2)   = 


b00+b102l1+b01Z21+bllZllz21+b212l2z21 

-1      -1      -1  -1      -2      -2  -1 
1  a10ai  "a0laz  "allZl  Z2  "a10Zl   a21Zl  Z2 


B(zl'z2} 
1-A(z1/Z  ) 


(IV. 1) 


Our  problem  will  be  drawing  a  detailed  signal  flowgraph  for 
the  system  function  H(z,,z_).   We  can  do  this  simply  enough  by 
combining  the  flowgraphs  on  Figures  2-2  and  2-3  to  get  the 
1NPOT     .  /T\  Poo    -f2\ i /X^ i _vX\    0^»«T 


Figure  4-1 

This  flowgraph  can  be  made  even  simpler  because  the  shift  opera- 
tion is  distributive  over  addition.   We  can  combine  the  two  z_ 
operators  into  a  single  one,  yielding  the  following  flowgraph. 
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IMP^T 


OUTPUT 


Figure  4-2 
Doing  so  reduces  the  number  of  shift  operators  that  need  to  be 
implemented  and  consequently  the  amount  of  storage  necessary. 
There  are  other  signal  flowgraphs  which  give  rise  to  the 
desired  system  function  H(z,,z  ).   For  example,  we  could 
invert  the  order  of  the  B-fz^z  )  filter  and  the  feedback 
loop  containing  A  (z,,z  )  to  obtain  the  block  diagram: 


\n  puT 


Afcti 


B&2) 


owtWT 


Figure  4"- 3 


Then,  when  we  substitute  Figures  2-2  and  2-3  for  the  blocks  as 
before,  the  two  z,   chains  will  contain  the  same  data  and  can 
be  merged  to  yield  the  signal  flowgraph  in  Figure  4-4.   This 
glowgraph  has  a  total  of  four  shift  operators,  and  it  mini- 


mizes the  number  of  z,       operators. 
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\N9WT 


QUT»<^T 


Figure  4-4 

Another  signal  flowgraph  that  minimizes  the  number  of  z, 
operators  may  be  obtained  from  Figure  4-4  by  the  2-D  transposition 
theorem  to  obtain  a  transposed  network.   Like  i-ts  1-D  counter- 
part [Ref.  9],  the  2-D  transposition  theorem  states  that  the 
transposed  network,  which  is  obtained  by  reversing  the 
directions  of  all  the  arrows  in  a  signal  flowgraph,  will  have 
the  same  system  function  as  the  original  network.   If  we 
reverse  the  direction  of  all  the  arrows  in  Figure  4-4  and  then 
redraw  the  flow  graph  with  the  input  port  on  the  left  and  the 
output  port  on  the  right,  we  get  the  flowgraph  shown  in  Figure  4-5 


OUTPUT 


Figure  4-5 
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This  transposed  flowgraph  may  be  preferred  in  implementations 
with  limited  wordlengths  since  the  attenuation  due  to  the 
"zeros"  of  H(z,,z  )  occurs  before  the  gain  due  to  the 
"poles"  thus  lessening  somewhat  the  possibility  of  arithmetic 
overflow  in  the  intermediate  computations. 

Using  the  notion  of  transposition  at  both  the  flowgraph 
level  and  the  block  diagram  level  (note  that  Figure  2-2  is  the 
transpose  of  Figure  2-1)  the  flowgraph  can  be  manipulated  to 
yield  a  realization  that  minimizes  the  total  number  of  shift 
operators . 

As  we  saw  earlier,  however,  a  z„   operator  will  require 
substantially  more  storage  than  a  z  1  operator  for  a  row-by- 
row  ordering  of  input  samples.   Consequently,  it  may  be  more 
economical  to  minimize  not  the  total  number  of  shift  operators 
(as  in  the  1-D  case)  but  the  number  of  z„   operators. 

If  the  filter  is  realized  by  using  a  separate  micro- 
processor to  compute  samples  of  each  node  signal,  storage 
may  be  less  of  an  issue. 

In  this  case,  we  may  want  to  minimize  the  total  number  of 
nodes  in  a  flowgraph  in  order  to  reduce  the  number  of  micro- 
processors in  an  implementation. 

As  digital  technology  progresses,  the  relative  costs  of 
storage,  computation,  and  interconnectivity  keep  changing.   In 
the  future  digital  systems  designers  may  have  radically  differ- 
ent criteria  for  optimizing  a  filter  realization. 
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B.   A  SECOND  ORDER  MODEL 

Looking   at  the  flowgraph  in  Figure  4-5  and  developing  a 
state  variable  implementation  from  it,  we  shall  call  the  output 

of  the  top  z,   operator  R,  (i,j)  ,  the  output  of  the  lower 

—  1  — 

z,   operator  R2  ( i , j  )  ,  the  output  of  the  left  z  1  operator 

S,(i,j)  and  the  output  of  the  right  z„   operator  S-(i,j)  as 

indicated: 


Sfiti 


S.Ci.J) 


Figure  4-6 


H(zrz2 


b00+b10Zl1+b01Z21+bllZllz21+b21Zl2z21 
,      -1      -1      -1  — 1      -2      -2  -1 
i"a10Zl  "a01Zl  "allZl  z2   a20Zl  "a21Zl  Z2 


R1(i+l,j) 
R2 (i+1, j) 

s,  (i> j+1) 

S2 (i,j+l) 


■10 


1  bii+boiaio  au+aoiaio 


a2Q   0   b21+bQ1al0   a21+a01a2Q 


01 


taU/J) 


R2 (i, j) 
S1(i,j) 


s2(i,j 


m 
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bio+booaio 

b00a20 


00 


u(i,: 


58 


(IV. 2) 


Y(i,j)   = 


[1   °   b0l   a01 


R2(i' J) 
S1(i,j) 


S2  Ci,j 


+   [bQ0]u(i,j)   (IV. 3) 


Defining 


bll  =  bll  +  aioboi 


all  "  an  +  aioaoi 


b21   =   b21  +  a20bQ1 


a2i   "   a21  +  a20a01 


In  general  the  foregoing  equations  can  be  written  as: 


'ij   =   bij+ai0b0i 


a.  .   =   a .  .+a . na.  . 
ij       i]   iO  0] 


Now  we  can  give  an  expanded  version  of  (IV- 2 
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Equations  (IV. 2)  and  (IV. 3)  represent  an  algorithm  for  comput- 
ing the  samples  of  the  output  signal  from  the  samples  of  the 
input  signal.   Just  as  in  the  preceding  subsection,  the  amount 
of  memory  required  to  store  the  state  variables  depends  on  the 
order  in  which  the  output  samples  are  to  be  computed.   It  is 
possible  to  envision  a  multiprocessor  architecture  for  comput- 
ing equation  (IV. 4)  by  assigning  each  processor  the  responsi- 
bility of  computing  the  next  value  of  a  particular  state 
variable  given  the  current  input  value  and  the  current  state- 
variable  values.   Equation  (IV. 3)  could  be  implemented  by  a 
filter  microprocessor  to  generate  the  desired  output  signal 
values . 

In  such  an  architecture,  minimization  of  the  number  of 
microprocessors  corresponds  to  the  minimization  of  the  number 
of  state  variables,  a  problem  studied  thoroughly  in  the 
literature.   Other  state-variable  forms  with  the  same  number 
of  state  variables  can  also  be  found  that  will  realize  the 
same  system  function  H(z, fz  )  and  may  exhibit  lower  coefficients 
of  sensitivity  or  round-off  noise  [Refs.  2,10]. 

For  the  special  case  of  "all-pole"  2-D  IIR  filters,  that  is, 
filters  with  a  system  function  of  the  form: 

H(z,  ,z0) 


A(z1,z2) 


where  bnQ  is  a  constant  and  AU-^z^  is  a  2-D  polynomial,  it 
can  be  shown  that  state  variable  realizations  based  on  signal 
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flowgraphs,  using  the  output  of  the  shift  operators  as  the 
state  variables,  require  the  minimum  number  of  state  variables. 
They  are  minimal  realizations  [Ref.  2]. 

From  the  above  equations  corresponding  to  the  second  order 
Roesser  model,  the  program  in  Appendix  C,  was  written.   This 
program  uses  the  values  of  coefficients  of  H(z,,z  )  as  inputs 
and  it  generates  an  output,  y(i,j).   Next,  the  program  finds 
the  2-D  Fourier  transform  of  this  output  matrix,  and  compares 
it  to  the  transfer  function  H  (co,  ,  oj.-,  )  . 
Numerical  Example 

In  the  following  three  examples  (first  and  second  orders), 
we  use  the  coefficients  of  first  and  second  order  transfer 
functions.   We  consider  the  special  case  of  "all-pole"  2-D  IIR 
filters,  i.e.,  filters  with  a  transfer  function  of  the  form: 


E(,„.,l 


A(z1,z2)      A(z1,z  ) 


where  bnn  is  constant  (unity  in  our  case)  and  A(z, ,z  )  is  a 
2-D  polynomial.   It  can  be  shown  that  state  variable  realiza- 
tions based  on  signal  flowgraphs,  using  the  output  of  the  shift 
operators  as  the  state  variables ,  require  the  minimum  number  of 
state  variables.   They  are  minimal  realizations  [Ref.  2]. 

For  the  third  program  we  have  a  graph  for  the  case  of  a 
BP  filter. 
Example  #4 

H(z,,z  )   =   — ^-- (IV. 6) 

1-0. 2zx  -0.5z2  -0.2z1  z2L 
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The  |Y(m,n) |  for  this  example  is  plotted  in  Fig.  4-la.  The 
corresponding  contour  map  of  2-D  surface  is  shown  in  Figure 
4-lb. 

Example  #5 

H(Z1'Z2)   =  -■> =7 =1T=I =T^T    (IV'7) 

l-0.Z5z11-0.345z21-0.125z1  z2  -0.1z1  z2 

The  2-D  D.F.T. |y(m,n)|  of  the  output  of  this  filter  is  shown 
in  Fig.  4-2a.   The  corresponding  contour  map  is  shown  in 
Fig.  4-2b. 
Example  #6 

-0.12  5+0.2  5z~1+0.125z~l-0.125z~1z~1+0.12  5z~2z~1 

H(z1,z2)   = -r— j 

1  +  z,  z? 

1  (IV. 8) 

|Y(m,n)|  for  this  examle  and  the  corresponding  contour  map  are 
shown  in  Figs.  4-3a  and  4-3b,  respectively. 

For  reasons  of  verification,  as  before,   Y(m,n)   was  com- 
pared to  the  actual  transfer  function   H(uk,uO   for  examples 
4,5,6.   These  transfer  functions  plots  and  the  corresponding 
contour  maps  are  shown  in  Fig.  4-4a,b,  Fig.  4-5a,b  and  Fig. 
4-6a,b  for  the  examples  4,5,6,  respectively. 

C.   EXTENSION  OF  THE  2-D  STATE  SPACE  MODELS  TO  HIGHER  ORDER 
TRANSFER  FUNCTIONS 

1 .   Introduction 

During  recent  years,  several  authors  (Attasi  [Ref.  11, 

Fozmasimi  and  Mazchesini  [Ref.  13],  Givone  and  Roesser  [Ref. 
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17] )  have  proposed  different  state  space  models  for  2-D 
systems.   They  have  also  suggested  some  extensions  of  the 
usual  1-D  notions  of  controllability,  observability,  and 
minimality  to  the  2-D  case. 

However,  these  results  are  not  quite  satisfactory. 
They  either  lack  motivation  for  the  state-space  models  intro- 
duced or  the  notion  of  state-space  is  improperly  defined.   In 
Chapter  II  we  started  with  a  comparison  of  all  the  current 
models  based  on  a  practical  (circuit-oriented)  point  of 
view  and  on  a  proper  definition  of  state.   It  is  shown  that 
the  model  of  Roesser  is  the  most  satisfactory,  in  that  it  is 
also  the  most  general  since  the  Attasi  and  Fozmasimi  Mazchesimi 
models  can  be  imbedded  in  the  Givone  and  Roesser  model. 

In  Chapter  II  we  pointed  out  that  a  major  difference 
between  1-D  and  2-D  systems  is  that  in  the  2-D  case  a  global 
state  (which  preserves  all  past  information)  and  a  local  state 
(which  gives  us  the  size  of  the  recursions  of  the  2-D  filter) 
can  be  introduced. 

2 .   Extension  for  2-D  Systems 

In  [Ref.  14],  Fozmasimi  and  Mazchesimi  use  the  algebraic 
point  of  view  of  "Nerode"  equivalence.   In  this  framework,  the 
state  space  arises  from  the  factorization  of  the  2-D  input/ 
output  map.   Fozmasimi  and  Mazchesini  were  the  first  to  realize 
that  a  major  difference  between  1-D  and  2-D  systems  is  that  we 
can  introduce  a  global  state  and  a  local  state  in  the  2-D  case. 

The  global  state  (which  is  of  infinite  dimensions,  in 
general)  preserves  all  the  past  information  while  the  local 


76 


state  gives  us  the  size  of  the  recursions  to  be  performed  at 
each  step  by  the  2-D  filter.   However,  Fozmasimi  and  Mazchesini 
failed  to  exploit  fully  the  structure  of  the  global  state  and 
its  relation  to  the  local  state,  so  that  the  state  space  model 
they  introduced  is  unsatisfactory  in  the  sense  that  what  they 
introduce  as  the  state  is  really  only  a  "partial  state"  (as 
defined  by  Wololich  [Ref .  15]  for  1-D  systems) .   Indeed,  this 
partial  state  does  not  obey  a  first-order  difference  equation 
(the  notion  of  first  order  difference  equation  for  linear  sys- 
tems or  partially  ordered  sets  has  been  defined  by  Mullans  and 
Elliot  in  [Ref.  16] \      Attasi's  model  suffers  from  the  same 
drawback  as  the  Fozmasimi  and  Mazchesini  one. 

On  the  other  hand,  Givone  and  Roesser  [Refs.  17,18,1] 
have  used  a  "circuit  approach"  to  the  problem  of  state  space 
realization  for  2-D  systems.   They  present  a  model  in  which 
the  local  state  is  divided  into  a  horizontal  and  a  vertical 
state  which  are  propagated,  respectively,  horizontally  and 
vertically  by  first-order  difference  equations.   From  this 
point  of  view  the  global  state  appears  as  the  boundary  condi- 
tion necessary  to  propagate  the  state-space  equations. 

However,  Roesser  did  not  provide  much  motivation  for 
the  introduction  of  such  a  model  and  seemed  unaware  of  the 
full  circuit  interpretation  of  their  model  since  they  were  not 
able  to  implement  an  arbitrary  2-D  transfer  function,  say 

b(z, , z  ) 

H(z,  ,zj 


'1'  2'  a(zx,z  ) 
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Mitra  et  al  gave  an  answer  in  [Ref.  19]  by  presenting 
an  implementation  method  for  2-D  transfer  functions  using  delay 
elements  z,-*-  and  z~  .   We  shall  see  below  that  this  approach  is 
consistent  with  Roesser's  model.   It  is  shown  in  [Ref.  8]  that 
Roesser's  model  appears  naturally  as  a  way  to  describe  the  local 
state  properties.   For  a  (n,m)  2-D  transfer  function, 


n    m 

■ L  ~      Ln      11  1 


"VD 


b(zl'Z2)      i=0  i=0   ij  X   2 
H(ZwZ0)   =    ;  x  *.       =      X  "  J" (IV. 9 


a(z,,z  )       n   m 

y     y    a.  .zn 

i=0  j-0   l:  x 


"V3 


exhibits  some  canonical  state-space  forms  (controllability, 
observability),  which  can  also   be  written  as, 


.1   ^(z^z;1 

H(zx,z2)   =   ±-5 (IV. 10) 

I    ai(z-I  )zT1 

i=0 


Without  loss  of  generality,  we  can  assume  an0  =  1  and  we 
denote 


a0(zt  )   =   1  +  a0(za  ) 


Thus,  using  1-D  realization  technique,  H(z,,z  )  of  Eq.  (IV. NO 
can  be  used  as  shown  below  in  Fig.  4-7. 
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Figure  4-7 


The  realization  is  almost  achieved:   in  addition  to  the 
n-horizontal  delay  elements,  we  need  only  m  vertical  delay 
elements  to  implement  the  feedback  gains  (a-(z_  ),  i  =  0,1,..., m) } 
and  m  other  vertical  delay  elements  to  implement  the  readout 
gains  {b.(z2  ),  i  =  0,1,..., m}.   Thus  the  complete  realization 
shown  in  Fig.  4-8  requires  only  n+2m  dynamic  elements.   This 
realization  is  a  standard  (canonical)  one;  its  structure  is  very 
simple  and  it  involves  only  real  gains.   Note  also  that  we 
need  fewer  dynamic  elements  than  was  suggested  by  the  imple- 
mentations of  [Ref .  19] . 

As  mentioned  in  Section  (b) ,  circuit  implementations 
with  delay  elements  z,   and  z~   are  in  a  one-to-one  corres- 
pondence with  state-space  models  of  Roesser's  type.   The 
outputs  of  the  z,1  delays  are  the  horizontal  states  and  the 
outputs  of  the  z~   delays  are  the  vertical  states.   Thus  the 
implementation  of  the  following  figure  can  be  transformed 
readily  into  the  following  state-space  model. 
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R(.i+l,j) 

S1(i,j+1) 
S2(i,j+1) 


A 


y(i,j)      =     c 


R(i,j) 

S1Ci,j) 

2, .     ., 
S     (i,D) 


R(i,j) 

S(i,j) 


+   bu(i,j) 


(IV. 11 


where : 


C      =       [b.n    ...    b    n         -b_n       0    ...    0      1      0       ...       0] 
10  nO  00 


[1   0  ...  0    aQl  ...  aQm   bQ1  ...  bQm]   (input  vector) 
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The   expanded    form  of    Eq.    IV-11    can   now   be    shown    as 
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D.   PROGRAM  AND  EXAMPLES  FOR  ROESSER'S  EQUATIONS  USING 
KUNG'S  MODEL 

This  program  (Appendix  D)  takes  as  initial  conditions  one 
horizontal  state  and  two  vertical  states.   The  order  of  hori- 
zontal states  is  given  by  N  and  the  order  of  the  vertical 
states  by  M. 

We  give  two  examples,  one  for  N  =  2  and  M  =  2  (Example  7) 
(two  orders  for  horizontal  states  and  2  orders  for  vertical 
states)  and  one  for  N  =  4  and  M  =  3  (Example  8)  (four  orders 
for  horizontal  states  and  three  orders  for  vertical  states) . 
The  first  example  is  for  a  matrix  2x2  and  the" second  example 
4x4. 
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E.   NUMERICAL  EXAMPLES  FOR  KUNG ' S  MODEL 

The  following  presents  three  examples.   The  first  one 
corresponds  to  an  "all-pole"  2-D  low-pass  filter.   The  second 
one  is  an  "all-zero"  2-D  band-pass  filter  (sin  w,  sin  w  ) . 
The  third  one  is  also  a  band-pass  filter.   All  these  examples 
are  second  order.   The  outputs  of  these  examples  are  produced 
using  Kung's  [Ref.  8]  state-space  model.   In  this  formulation, 
for  a  second  order  system,  we  require  two  horizontal  states-- 
Rl(i,j)  and  R2(i,j)  and  four  vertical  states,  Sl(l)(i,j), 
Sl(2)(i,j),  S2(l)(i,j)  and  S2(2)(i,j).   The  program  listing 
for  implementing  this  model  is  given  in  Appendix  D. 
Example  #9 

The  system  parameters  and  the  initial  conditions  chosen 
for  this  example  are  as  listed  in  Table  4.1.   The  2-D  D.F.T. 
|Y(m,n) |  of  the  output  sequence  y(i,j)  produced  by  the  program 
in  Appendix  D  is  shown  in  Fig.  4-9a.   The  corresponding  contour 
map  is  shown  in  Fig.  4-9b. 
Example  #10 

The  parameter  coefficients  and  the  initial  conditions  for 
this  example  are  listed  in  Table  4.2.   The  2-D  D.F.T. 
sequence  |Y(m,n) |  for  this  example  is  illustrated  in  Fig.  4-10a, 
and  Fig. M.  10b  shows  the  associated  contour  map. 
Example  #11 

The  parameter  coefficients  and  the  initial  conditions  for 
this  example  are  listed  in  Table  4.3.   The  2-D  D.F.T.  sequence 
|Y(m,n) |  for  this  example  are  illustrated  in  Fig.  4-lla  and 
Figure  4-llb  shows  the  associated  contour  map. 
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TABLE  4.1 


NUMBER  OF  HQRIZONTGL  STATES ( N=lto4> :  2 


NUMBER  OF  VERTICAL  STATES <M= 1 to4) :  3 


DIMENSION  OF  OUTPUT ( It o25) 


15 


ENTER  INITIAL 

CONDITIONS 

R  1  (  1, 

1)  : 

0 

R  3(1, 

1)  : 

o 

R  1(1, 

2)  : 

0 

R  3(1, 

2)  : 

0 

R  1(1, 

3)  : 

0 

R  2(1, 

3)  : 

0 

R  1  (1, 

4)  : 

0 

R  2(1, 

4)  : 

0 

R  1(1, 

5)  : 

0 

R  2<  I, 

5)  : 

0 

R  1(1, 

6)  ■ 

0 

R  2(1, 

S)  : 

0 

R  1  (  1, 

7) 

0 

R  2(1, 

7)  ■ 

0 

R  1  (  1, 

8) 

0 

R  2(1, 

a> 

0 

R  1  (1, 

9) 

•  0 

R  2(1, 

S) 

o 

R  1  (1, 

10) 

:  0 

R  2(1, 

10) 

:  0 

R  1  (  1, 

id 

:    0 

R  3(1, 

1 1) 

:  0 

R  1  (  1, 

12) 

:  0 

R  3(1, 

12) 

:  0 

R  1  (  1, 

13) 

:  0 

R  3<  1, 

13) 

:  0 

R  1  (  1, 

14) 

:  0 

R  2(1, 

14) 

:  0 

R  1  (  1, 

15) 

:  0 

R  3(1, 

15) 

:  0 

ENTER 

INI 

riAL 

CONDITIONS 

SI  (  1) 

(  1 

.  1)  : 

0 

SI  (  2) 

(  1 

,  1)  : 

0 

SI  (  1) 

(  2 

,  1)  : 

0 

SI  (  3) 

(  2 

,  1)  : 

0 

SI  (  1) 

(  3 

,  1)  : 

0 

SI  (  3) 

(  3 

,  1)  : 

0 

SI  (  1) 

(  4 

,  1)  : 

0 

SI  (  2) 

(  4 

,  1)  : 

0 

SI  (  1) 

(  5 

,  1)  : 

o 

SI  (  2) 

(  5 

,  1)  : 

0 

SI  (  1) 

(  S 

,  1)  : 

0 

SI  (  2) 

(  S 

,  1)  : 

0 

SI  (  1) 

(  7 

,  1)  : 

0 

SI  (  2) 

(  7 

,  1)  : 

0 

SI  (  1) 

(  a 

,  1)  : 

0 

SI  (  2) 

(  a 

,  1)  : 

0 

SI  (  1) 

(  3 

,  1)  : 

0 

SI  (  2) 

(  s 

,  1)  : 

0 

SI  (  1) 

(10 

,  1)  : 

0 

SI  (  2) 

(10 

,  1)  : 

0 

SI  (  1) 

(ii 

,  1)  : 

0 

SI  (  2) 

(ii 

,  1)  : 

0 

SI  i     1) 

(12 

,  1 )  : 

0 

FOR  HORIZONTAL  R(#.*) 


FOR  VERTICAL  SI (#.*) 


31  ( 

I)   ( 12,  l>  : 

0 

SI  ( 

2)  ( 13,  1)  : 

0 

31  ( 

1)  (  W,  l  )  : 

o 

SI  ( 

2)  ( 14,  l )  : 

0 

SI  ( 

1)  ( 15,  1)  : 

o 

SI  C 

2)  (15,  1)  : 

0 

ENTER  INITIfiL 

CONDITIONS 

S2< 

1)  (  .1,  1)  : 

0 

32  < 

2)  (  1,1): 

0 

S2< 

1)  (  2,  1)  : 

0 

S2< 

2)  (  2,  1)  : 

0 

S2( 

1)  (  3,  1)  : 

0 

32  ( 

2)  (  3,  1)  : 

0 

S2( 

1)  (  4,  1)  : 

0 

S2< 

2)  (  4,  1)  : 

0 

32  ( 

1)  (  5,  1)  : 

0 

S2< 

2)  (  5,  1)  : 

0 

S2< 

1)  (  S,  1)  : 

0 

S2< 

2)  (  S,  1)  : 

0 

S2( 

1)  (  7,  1)  : 

0 

S£( 

2)  (  7,  1)  : 

0 

32  ( 

1)  (  3,  1)  : 

0 

S2( 

2)  (  8,  1)  : 

0 

S2( 

1)  (  3,  1)  : 

0 

32  C 

2)  (  9,  1)  : 

0 

S2( 

1)  (10,  1)  : 

0 

32  ( 

2)  (10,  1)  : 

0 

S2( 

1)  (  11,  1)  : 

0 

S2t 

2)  (11,  1)  : 

0 

S2( 

1)  (12,  1)  : 

0 

32  ( 

2)  (12,  1)  : 

0 

S2( 

1)  (12,  1)  : 

0 

S£( 

2)  (13,  1)  : 

0 

S2( 

1)  (  14,  1)  : 

0 

S2C 

2)  (14,  1)  : 

0 

S2( 

1)  (15,  1)  : 

0 

32  ( 

2)  (15,  1)  : 

0 

FOR  VERTICPL  S2(#.#) 


ENTER  VPLUES  FOR  THE  INPUT  VECTOR (#.  *) 
a<0  1) :  -O. 35 
a(0  2) :  0 
b  ( 0  1 )  :  0 
b(0  2)  :  0 

3  OF  THE  TRPNSITION  MBTRI X ( #. #> 


ENTER  EL 

EMENT3 

a< 

10)  : 

-0. 125 

a( 

20)  : 

-0.  25 

a( 

1 

1)  : 

-0.  1 

a( 

•3 

1)  : 

0 

a( 

1 

2)  : 

0 

a( 

2 

2)  : 

-0.  1 

b( 

1 

1)  : 

0 

b( 

o 

1)  : 

0 

b( 

1 

2)  : 

0 

b( 

o 

2)  : 

0 

ENTER  VALUES  FOR  THE  OUTPUT  VECTOR  <**.#> 

b(00) :  1 
b(  10)  :  0 
b(  20) :  O 
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TABLE  4.2 

NUMBER  OF  HORIZONTAL  STATES (N= It o4> :  2 
NUMBER  OF  VERTICAL  STATES (M=lto4) :  2 

DIMENSION  OF  OUTPUT ( lto25) :  17 


ENTER  INITIAL 

CONDITIONS 

R  1  ( 1,   1)  : 

Q 

R  2(1,  1)  : 

0 

R  1(1,  2)  : 

0 

R  2(1,  2)  : 

0 

R  1(1,  3) : 

0 

R  2(1,  2)  : 

0 

R  1(1,  4)  : 

0 

R  2(1,  4)  ; 

0 

R  1(1,  5) : 

0 

R  2(1,  5)  : 

0 

R  1(1,  S>  • 

0 

R  2(1,  6)  • 

0 

R  1(1,  7) 

o 

R  2(1,  7) 

0 

R  1(1,  3) 

0 

R  2 (  1,  8) 

0 

R  1(1,  9) 

0 

R  2(1,  9) 

0 

R  1  ( 1,  10) 

.  0 

R  2(1, 10) 

0 

R  1(1,11) 

:  0 

R  2(1,11) 

:  0 

R  1 (1, 12) 

:  0 

R  2(1, 12) 

:  0 

R  1 (1, 13) 

:  0 

R  2(1, 13) 

:    0 

R  1 ( 1, 14) 

:  0 

R  2T1, 14) 

:  0 

R  1 (1,  15) 

:  0 

R  2(1, 15) 

:  0 

R  1 ( 1, 16) 

:  0 

R  2(1, IS) 

:  0 

R  1 (1, 17) 

:  0 

R  2(1,17) 

:  0 

ENTER  INT 

flAL  CONDITIONS 

SI  (  1)  (  1 

,  1)  : 

0 

SI  (  2)  (  1 

,  1)  : 

0 

SI (  1)  (  2 

,  1)  : 

0 

SI  (  2)  (  2 

.  1) 

0 

SI (  1)  (  3 

,  1)  : 

0 

SI  (  2)  (  3 

,  1) 

0 

SI (  1)  (  4 

,  1> 

0 

SI  (  2)  (  4 

,  1) 

0 

SI (  1)  (  5 

,  1) 

0 

SI (  2) (  5 

,  1> 

0 

SI (  1)  (  6 

,  1) 

0 

SI  (  2)  (  6 

,  1) 

:  0 

SI (  1)  (  7 

,  1) 

!  0 

SI (  2) (  7 

,  1) 

:  0 

SI (  1)  (  a 

,  1) 

■  0 

si (  2)  (  a 

,  1) 

:  0 

Si (  1)  (  9 

,  1) 

:  0 

SI (  2) (  9 

,  1) 

:  0 

S 1  (  1 )  ( 1 0 

,  1) 

:  0 

3 1  (  2 )  (  1  0 

,  1) 

:  0 

SI (  1) (11 

,  1  ) 

:  0 

FOR  HORIZONTAL  R<#.#> 


FOR  VERTICAL  SI (#. #) 
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31  ( 

2) (12,  U 

0 

Si  ( 

1)  (13,  1)  . 

0 

SI  ( 

3)  (13,  1) 

0 

SI  ( 

1  )  (  14,  1) 

0 

31  ( 

2)  (14,  1) 

0 

SI  ( 

1)  (IS,  1) 

0 

SI  ( 

2)  (IS,  1) 

0 

31  ( 

1)  ( IS,  1) 

0 

SI  ( 

2)  (16,  1) 

:  0 

31  ( 

1)  (17,  1) 

0 

31  ( 

2) (17, 1) 

:  0 

ENTER  INITIAL  CONDITIONS  FOR  VERTICAL 

S£< 

1)  (  1,  1) 

:  0 

S2( 

3)  (  1,  1) 

:  0 

32  < 

1)  (  2,  1) 

:  0 

S£( 

2)  (  2,  1) 

:  0 

S£< 

1)  (  3,  1) 

:  0 

S£< 

2)  (  3,  1) 

:  0 

S2( 

1)  (  4,  1) 

:  0 

S£< 

3)  (  4,  1) 

:  0 

S£< 

1)  (  5,  1) 

:  0 

S£( 

2)  (  5,  1) 

:  0 

S£( 

1)  (  S,  1) 

:  0 

S£< 

2)  (  S,  1) 

:  0 

S£< 

1)  (  7,  l) 

:  0 

S£< 

2)  (  7,  1) 

:  0 

S2< 

1)  (  S,  1) 

:  0 

S£( 

2)  (  a,  i) 

:  0 

S£( 

1)  (  9,  1) 

:  0 

S£< 

3)  (  9,  1) 

:  0 

S£( 

1)  (10,  1) 

:  0 

S£< 

2)  (10,  1) 

:  0 

S£< 

1)  (11,  1) 

:  0 

S£< 

2)  (11,1) 

:  0 

S2C 

1) (12, 1) 

:  0 

S2  '. 

2)  (12,  1) 

:  0 

S£< 

1)  ( 13,  1) 

:  0 

S2( 

3) (13, 1) 

:  0 

S£( 

1) ( 14, 1) 

:  0 

S£( 

2) (14, 1) 

:  0 

S2( 

1) (15, 1) 

:  0 

S£< 

2)  (15,  1) 

:  0 

S2< 

1)  (16,  1) 

:  0 

S£< 

2)  (16,  1) 

:  0 

S£< 

1) (17, 1) 

:  0 

S£< 

2) (17, 1) 

:  0 

;2(*.  *) 


ENTER  VALUES  FOR  THE  INPUT  VECTOR (#. #> 
a  ( O  1 )  :  O 
a(0  2) :  0 
b  ( 0  1 )  :  O 
b(0  2) :  0. 125 


ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX (#. #) 

a(  10)  :  0 

a (  20 )  :  0 

a(  1  1)  :  0 

a  (  2  1 )  :  0 

a(  1  2) :  0 

a(  2    2):     0 

b(  1  1)  :  0 

b  (  3  1  )  :  0 

b<  1  3) :  0 

b(  2    2):     -O. 

125 
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ENTER  VALUES  FOR  THE  OUTPUT  VECTOR  ( •*.  #) 

b (00) :  0. 133 
b(  10)  :  O 
b(  £0) :  0. 135 


*****  INPUT  VECTOR 
. .  00       . 00       . 00       . 00       .  00 


♦***♦  OUTPUT  VECTOR  ***** 
00      .12     -.13      .00     1.00      .00 

*****  TRANSITION  MATRIX  ***** 


.  00 

.  00 

-1.  00 

.00 

.  00 

.  00 

1 .  00 

.  00 

.  00 

.  00 

.  00 

.00 

.  00 

.  00 

.  00 

1 .  00 

.  00 

.  00 

.00 

.  00 

.  00 

.00 

.  00 

.00 

.00 

.  00 

.  00 

.  00 

.  00 

1.  00 

.00 

-.  13 

-.  13 

.  00 

.00 

.00 
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Once  again,  to  verify  the  correctness  of  our  program, 
the  D.F.T.  |Y(m,n)|  was  compared  to  |  H  (go,  ,  u)-)  |  .    H(a),  ,wJ 
and  the  corresponding  contour  maps  are  shown  in  Fig.  4-12a,b, 
Fig.  4-13a,b  and  Fig.  4-14a,b  for  examples  9,  10  and  11, 
respectively. 

F.   SUMMARY  OF  PROGRAMS  DEVELOPED 

The  programs  which  have  been  written,  cover  the  following 
orders  based  upon  the  different  models. 

Appendix      Order  Model  #  of  States 

A  1st  Roesser  1  horizontal,  1  vertical 

C  2nd  Roesser  2  horizontal,  1  vertical 

D  Multi-order  Kung  v\  horizontal,  2vi  vertical 

In  order  to  check  the  program  listing,  the  same  first 
order  example  was  used  on  all  programs.   Identical  results 
were  obtained.   Similarly,  identical  second  order  examples 
were  used  in  Programs  C  and  D  and  produced  identical  outputs. 
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V.   USE  OF  SSPACK  PACKAGE 

A.   SSPACK 

"SSPACK"  is  a  "state  space  system  package,"  [Ref.  21]  that 
is  an  interactive,  state-of-the-art,  software  package  for  the 
analysis,  design,  and  display  of  one-dimensional  state-space 
systems.   The  work  which  follows  adapts  this  program  so  that 
it  can  be  used  to  produce  2-D  data  fields  from  state  space 
formulations.   A  brief  description  of  SSPACK  follows. 

SSPACK  is  useful  for  a  variety  of  applications  in  signal 
processing  and  control  [Ref.  22].   The  package  consists  of  a 
supervisor  which  controls  the  operation  of  the  software  and  a 
set  of  independent  programs  which  communicate  using  disk  files. 
The  core  of  the  package  are  the  pre-  and  post-processors.   The 
state-space  pre-processor  (SSPREP)  program  aids  in  preparing 
files  for  the  individual  algorithm  programs.   [Refs.  23,24]   The 
state  space  post-processor  (SSPOST)  program  displays  and  analyzes 
the  output  from  the  algorithms.   SSPREP  prompts  with  a  series 
of  questions  in  a  menu  format. 

SSPOST  is  an  interactive  command-drive  processor.   It  is 
designed  to  help  interpret  the  output  of  the  various  SSPACK 
algorithms,  and  display  time  histories: 

A  is  the  Nx  by  Nx  state  transition  matrix; 
B  is  the  Nx  by  Nu  input  transition  matrix; 
C  is  the  Nz  by  Nx  measurement  matrix; 
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D  is  the  Nz  by  Nu  feedthru  matrix; 

W  is  the  Nx  by  Nw  process  noise  matrix; 

V  is  the  Nz  by  Nv  measurement  noise  matrix. 
The  SSPACK  works  in  multi-order  form,  using  the  transfer 
function  of  the  1-D  digital  filter. 


U 


Figure  5-1 


The  present  objective  is  to  use  SSPACK  with  a  2-D  input  data 
field  and  through  the  same  transfer  function,  1-D  digital 
filter,  to  accomplish  2-D  output  data  field. 

B.   DESIGN  OF  2-D  DIGITAL  FILTERS  USING  1-D  DIGITAL  FILTER 
STRUCTURES 

The  idea  of  using  two  types  of  dynamic  elements  is  not  very 
abstract;  it  is  very  natural  in  delay-differential  systems. 
However,  before  considering  its  practical  applications  to 
image  systems,  two  remarks  have  to  be  made.   The  first  is  be- 
cause the  "spatial"  dynamic  elements  seem  un implement able ,  and 
we  need  to  replace  them  by  time-delay  elements.   Secondly,  in 
order  to  have  a  finite  order,  we  shall  only  consider  a  bounded 
frame  system,  i.e.,  we  assume  that  the  picture  frame  of  interest 
is  an  M  x N  frame  (with  vertical  width  M  and  horizontal  length 
N) .   Note  that  in  order  to  use  time  delay  elements,  we  need 
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first  to  find  a  way  to  code   a  2-D  spatial  system  into  a  1-D 
(discrete  time)  system  and  vice  versa. 

Thus  we  propose  the  following  system,  composed  of  three 
subsystems  in  series: 

i)   The  Input  Scan  Generator  codes  the  2-D  spatial  input 
into  1-D  time  data  according  to  the  mapping  function 


t(-,*) 


t(i,j)   =   iM  +  jN  ,   0  <_   i    N-l 

0  <  j  <  M-l 


(V.l) 


where  M  and  N  are  relatively  prime  integers.   For  example,  we 
consider  a  2-D  input  data  u(i,j): 


u(i, j) 


(0,0) 
(1,0) 


(o,i: 
(i,i 


0,2)    (0,3 


(N-1,0)   (N-1,1) 


(0,M-1) 
(1,M-D 

(N-l, M-l) 


Scanning 

The  data  field  u(i,j)  is  scanned  to  produce  u(^)  as  follows 


u(i,  j ) 


(0,0),  (0,1),  (0,2),...,  (0,M-1)  ,(1,0),  (1,1), 
(1,M-1) , (N-1,0) . . . (N-l, M-l) 


(u(t) },  t  =  0,  1,  2,  M,  M+l,  M+2,  . . . ,  (M-l) (N-l) ,   t  =  iM  +  jN 


For  example, 
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u(i,j) 


10         0 
0         0         0 

0         0         0 


yields 


y(t)       =       [10      0      0      0...       0000 


ii)   A  1-D  (discrete  time)  digital  filter  processes  the 
1-D  data  generated.   This  subsystem  is  implemented  by  replacing 
z,   by  <5,  z~   by  A  in  a  2-D  circuit  realization  (e.g.,  2-D 
controller  form).   6  and  A  are  chosen  as: 


M 


6   =   D    =   M-units  delay  element 

N 
A   =   D    =   N-units  delay  element 


iii)   The  Output  Frame  Generator  decodes  the  1-D  (discrete 
time)  output  of  the  1-D  digital  filter  described  above  into 
a  2-D  (discrete-spatial)  picture  according  to  the  inverse 
mapping  of  (V.l) . 


(i(t) ,j (t) )   =   Pt  Mod  N,[t-(Pt  Mod  N)M]/N) 


(V.2 


where  P  is  a  unique  integer  such  that  PM-PN  =  1  and  0  <  P  <  N. 
This  formula  is  given  in  [Ref.  2].   Alternately,  we  can  compute 
(i, j)  as 

i   =   t  Mod  N 


and 


j   =   Quotient  (t/N! 
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For  example  we  suppose  t  =  19  with  N  =  10  and  M  =  9. 

19 
The  corresponding  value  in  the  2-D  case  will  be  i  =  Remainder {-=-q} 

19 
=  9  and  j    =   Quotient  {-=-q-}  =  1.   So  in  the  2-D  case  we  will 

have  (i, j)  =  (9,1) . 

Another  Example:   For  M  =  4  and  N  =  5,  the  single  index  i 

will  be  mapped  into  (i,j)  as: 


: 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

The  procedure  for  implementing  2-D  filters  using  1-D  filter 
structures  is  as  shown  below  in  Fig.  5-2. 


|N»VJT 


"Ci,i) 


F\£J_b 


I&ewfc- 


YQ-J) — *D.F.^ 


Figure  5-2 
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The  index  scanning  is  required  for  the  input  data  so  SSPACK 
can  be  carried  out  simply  because  the  input  is  assumed  to  be 
a  2-D  unit  pulse.   This  is  followed  by  implementing  the 
corresponding  1-D  filter  of  Fig.  5-3  using  SSPACK  to  convert 
the  1-D  output  from  SSPACK  to  a  program  for  output  index 
mapping--written  as  shown  in  Appendix  E.   The  2-D  Fourier 
transform  of  the  resulting  2-D  field  is  then  computed. 

Considering  a  bounded  frame  (M  xN)  system  it  is  interesting 
to  know  the  dimension  of  the  global  state  (or  initial  condi- 
tions) needed  to  process  the  M  xn  future  data  field.   Since 
vertical  states  convey  information  vertically,  all  the  verti- 
cal states  along  the  X-axis  are  necessary  initial  conditions 
and  their  dimension  is  mN.   Similarly,  all  the  horizontal  states 
along  the  Y-axis  are  necessary  initial  conditions  (with 
dimension  nM) .   They  convey  information  horizontally. 

Therefore,  in  the  bounded  frame  case  a  total  number  of 
mN+nM  are  needed  to  summarize  the  "past"  information.   This 
very  same  idea  can  be  used  again  from  a  computational  point 
of  view.   Indeed,  the  number  of  required  storage  elements  for 
recursive  computations  is  also  equal  to  mN+mN  if  initial  condi- 
tions are  not  zero.   However,  it  is  quite  often  the  case  that 
the  system  starts  with  zero  initial  conditions;  the  size  of 
storage  required  is  reduced  to  mN  (respectively,  nM)  which  is 
used  to  store  the  updated  data  row  by  row  (respectively,  column 
by  column) .   No  storage  is  needed  for  the  rest  of  the  initial 
conditions--nM  horizontal  states  (respectively,  mN  vertical 
states)  since  they  are  assumed  to  be  zero.   This  is  consistent 
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with  the  results  of  Read  [Ref.  24]  derived  from  a  direct 
polynomial  approach. 

Another  interesting  observation  concerns  the  dimension 
of  the  1-D  figital  filter  contained  by  our  2-D  digital  filter 
design  discussed  above.   Since  it  needs  nM  unit-delays  and 
mN-unit  delays,  the  corresponding  1-D  state-space  also  has  a 
dimension  equal  to  nM+mN.   Note  that,  despite  the  high  dimen- 
sion of  the  corresponding  1-D  filter,  its  high  sparsity  is  very 
encouraging  for  further  studies.   In  short,  following  the 
above  method  of  designing  a  2-D  filter,  for  the  first  order 
case, 


H<zrz2)   = _ -r-J.  (V.3 

1+a10Zl  +a01Z2  +allZl  Z2 


Using  the  above  approach  we  get  the  1-D  filter  realization  for 
this  2-D  filter  which  turns  out  to  be  as  shown  in  Fig.  5-3. 

The  detailed  matrix  equations  for  realizing  Eq.  (V.3) 
using  SSPACK  can  be  written  as,  The  SSPACK  produces  a  1-D 
sequence,  which  converted  into  a  2-D  sequence  using  the  output 
index  mapping  formulae  discussed  earlier.   The  listing  of  a 
program  which  does  this  mapping  is  shown  in  Appendix  E. 

After  obtaining  the  valid  2-D  output  data  sequence  y(i,j) 
we  next  compute  its  2-D  D.F.T.  to  produce  |Y(m,n) |  which  for  this 
example  is  plotted  in  Fig.  5-4a.   The  corresponding  contour 
map  is  as  shown  in  Fig.  5-4b. 
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VI.   CONCLUSIONS 

This  thesis  has  dealt  with  the  problem  of  modelling  2-D 
data  fields  in  the  state-space  domain.   First  of  all  we  have 
pointed  out  the  main  problems  associated  with  the  extension  of 
1-D  time-discrete  state-space  models  to  2-D  data  fields.   The 
remaining  part  of  the  thesis  has  been  divided  primarily  in  3 
parts . 

In  the  first  part  we  describe  Roesser's  [Ref.  5]  approach  to 
modelling  2-D  systems  in  the  state  space  domain.   Extensive  com- 
puter simulation  results  are  presented  to  verify  the  functioning 
of  this  approach.   This  modelling  approach  has  been  tried  out 
for  the  scalar  (1  x 1)  as  well  as  for  higher  order  (2  x 2)  etc., 
2-D  systems. 

The  second  part  deals  with  a  modification  of  Roesser's  approach 
as  described  by  Kung  [Ref.  7].   The  main  advantage  of  this  approach 
is  that  the  2-D  state-space  model  can  be  realized  as  a  2-D 
circuit.   More  importantly,  this  2-D  circuit  realization  can  be 
implemented  as  a  1-D  digital  filter.   Computer  simultation  studies 
that  have  been  carried  out  substantiate  the  making  of  this  model. 
The  1-D  filter  realization  obtained  in  this  part  turns  out  to 
be  a  very  convenient  starting  point  for  the  nezt  part  of  our 
effort,  dealing  with  the  use  of  the  1-D  SSPACK  commercial  soft- 
ware package  designed  for  dynamic  system  simulation. 

In  the  final  part  of  the  thesis,  we  make  use  of  the  1-D 
filter  realization  of  2-D  state-space  model  obtained  in  the  second 
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part,  and  implement  this  filter  using  SSPACK.   Some  additional 
programming  effort  reuiqred  for  input  and  output  mapping  was 
necessary.   Programs  for  converting  2-D  input  and  output  se- 
quences to  1-D  have  been  written  separately.   In  this 
fashion  we  have  succeeded  in  extending  the  applicability  of 
the  SSPACK  to  simulating  2-D  linear  systems  as  well.   Once 
again,  detailed  computer  simulations  have  been  carried  out 
to  verify  the  functioning  of  this  modification  of  the  SSPACK. 
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APPENDIX    A 


Line*     1 

1  3STCRAGH 

2  SLARGE 


naqe         I 
09-25-85 

12:3^:27 
Microsoft  F0RTRPN77  V2.  20  02/8** 


5 

7 

a 

3 
10 

11 

12 

12 
14 
15 
IS 
17 
18 
13 
20 
21 


*  * 

*  THE  PURPOSE  GF  THIS  PROGRAM  IS   TO   COMPUTE   PND  GRAPH  THE  * 

*  EQUATIONS  OF  ROBERT  P.  ROESSER  IN  THE  "DISCRETE  STATE-SPACE  * 

*  MODEL  FOR  LINEAR  IMAGE  PROCESSING".  * 

*  * 

*  EVANGELOS  THEOFILOU  * 

PROGRAM  2D-DATA-FIELD 

*-*♦-**  VARIABLE  DECLARATIONS  ***** 

REAL         R(25,  25) , S(25, 25)  ,  R 1 ( 2) , R2 (2) ,  Z (31 ,  3 1 ) , 
►  RLPART, IMGPART, ZF(31, 31) , VERTEX (15) , ZLEV(3l) 

INTEGER      MASK (3000) , LDIG (31) , LWGT (31) 
CHARACTER*!  ANSWER 
CHARACTER*20  CTEXT 


DATA 


XLOL/0. 0/, VLOL/0. 0/, XUPR/8. 5/, YUPR/7. 0/, 
ZLOUI/1.  0E25/,  IPROJ/0/,  NRNG/  100/ 


24  C 

25 

26  C 

27 

23 

29 

20 

31 

32 

34 

26 
37 

3a 

33 
40 

41 
42 
43 
44 
45 

46 
47 
48 
43 
50 
51 
52 
53 
54 

56 

57 
58 

53 


MAIN 


PROGRAM 


'HE  REuUIRED  VALUES  FOR  THE  MODEL  ***** 

HE  FOLLOWING  VAR I ABLES (*.*,.. ) 


INTER 
'  Al  : 


VALUES  -OR 


10  WRITE  (*. *)  ' 

WRITE  <*, 333) 

READ  (*, *)  Al 

WRI'E  (*,333)  'A2: 

READ  (*, *)  A2 

WRITE  (*, 333)  'A3: 

READ  (*, +)  A3 

WRITE  (*, 333)  ' A4: 

READ  (*, *)  A4 

WRITE  (*, 333)  '31: 

READ  (*, *)  Bl 

WRITE  (*, 333)  '32: 

READ  (*, *)  BE 

WRITE  (*, 333)  'CI: 

READ  (*, *)  CI 

WRITE  (*, 333)  'C2: 

READ  (*, *)  C2 
5   WRITE  (*, 402) 

READ   (*, *)  N 

IF  (N  . GT.  25)  GOTO  5 

WRITE  (*,  211)  'ENTER  '  ,  N,  '   INITIAL  CONDITIONS  FOR  MATRIX  R<#.  **)' 
DO  33  I  =  1 ,  N 

WRITE  (*,403)  'R(1,',I,'):  ' 
READ   (*, *)  R(l, I) 
99  CONTINUE 

WRITE  (*,211)  'ENTER  ' , N, '  INITIAL  CONDITIONS  FOR  MATRIX  S (#.#)' 
DO  100  I  =  1,N 

WRITE  (*,404)  'S(,',I,'D:  ' 
READ   (*. *)  S( I, 1 ) 
100  CONTINUE 
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Line# 

SO 

61 

62 

62 

64 

65  i 

66 

67 

66 

63 

70 

71 

7S 

73 

74 

75 

76 

77 

76 

73 

80 

31 

82 

83 

84 

85 

86 

87 

88 

83 

30 

31 

32 

33 

34 

35 

36 

37 

38 

33 

100 

101 

102 

103 

104 

105 

106 

107 

108 

103 

110 

111 

112 

113 

114 

115 

116 

1  17 

1  18 


Paqe    £ 
03-25-85 

12:34:2: 
7  Microsoft  FQRTRAN77  V3. 20  02/6- 

WRITE  (*, 413) 
READ   (*, 200)  ANSWER 
IF   ((ANSWER  .EQ.  '  Y» )  .OR.   (ANSWER  . EQ.  ' y'  )  )  GOTO  10 

U  =  1.0 

*****  COMPUTE  R  AND  S  MATRICES  ***** 

DO  101  I  =  1 , N 

DO  101  J  =  1 , N 

IF   (1+1  . LE.  N)   THEN 

R(I  +  1,J)  =  A1*R(I,J)  +  A2*S(I,J)  +  B1*U 
ENDIF 
IF   (J+l  . LE.  N)   THEN 

S(I,J+1)  =  A3*R(I,J)  +  A4*S(I,J)  +•  B2*U 
ENDIF 
U  =  0.  O 

101  CONTINUE 

****  FILL  0' s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 
DO  102  I  =  1,31 
DO  102  J  =  1,31 
Z(I,J)  =0.0 

102  CONTINUE 

*****  COMPUTE  Z  MATRIX  ***** 
DO  103  I  =  1,N 
DO  103  J  =  1, N 

Z(I,J)  =  Cl*R(J,I)  +  C£*S(J,I) 

103  CONTINUE 

*****  OUTPUT  THE  Z  MATRIX  ***** 

WRITE     (*,205)     '*****       Z       MATRIX       ' , N, '     X    ' , N, '        *****' 

WRITE  (*,212) 

DO  104  I  =  1,N 

WRITE  (*,  300)   (Z(I,J),  J  =  1,N) 

WRITE  (*,210) 

104  CONTINUE 
WRITE  (*,213) 

WRITE  (*, 418) 

READ   (*,200)  ANSWER 

IF   ((ANSWER  . NE.  ' Y' )  .AND.   (ANSWER  . NE.  ' y' ) )  GOTO  18 

*****  PSK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 
15   WRITE  (*, 210) 


WRITE  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 


*)  i  **♦   ENTER   PLOT   PARAMETERS  *** 

405) 

*)   AZIM 

406) 

*)   ELEV 

408) 

*)    ITRIM 

403) 

*)    IDIV 

411) 

133)   CTEXT 

401) 

200)   ANSWER 
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124 
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133 
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141 
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IF   <  (ANSWER  .ED.  'YM  .OR.   (ANSWER  .  EQ.  '  y'  )  )   THEN 

COLL  PLOTS (0,0, 2) 
ELSE 

CPLL  PLOTS ( 0, 33, 33 ) 
ENDIF 

CALL  WINDOW(XLOL, YLOL,  XUPR,  YUPR) 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS(Z,  31,31,N,N,AZIM,  EL2V,  0.  5,  0.  5,  7.  5,  5.  5,  ID IV,  0, 
*  3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (5. 5, 0. 3, 0. 2, ' AZIMUTH:     ',0.0,10) 

CALL  NUMBER (393. O, 333. 0, 0. 2, AZIM, 0. 0, 2) 

CALL  SYMBOL (5. 5, 0. 0, 0. 2, ' ELEVATION: ' , 0. 0, 10) 

CALL  NUMBER (333. O, 333. O, 0. 2, ELEV, 0. 0, 2) 

DY  =  (Z  (  1,  1) /30.  0)  *  ELEV 

CALL  P3D2D ( 1. 0,  1. 0,  Z  < 1,  1) -DY,  XR, YR) 

CALL  SYMBOL (XR, YR, 0. 25, ' *' , O. 0, 1) 

CALL  SYMBOL ( 1.  0, 0.  1, 0. 2,  ' *  =  ORIGIN'  , 0. 0,  10 ) 

CALL  SYMBOL (1. 0, 6. 75, O. 25, CTEXT,  0. 0,  20) 

CALL  SYMBOL (6. 0, 6. 5, O. 2, ' 2-D  DATA  FIELD' , 0. 0, 14 ) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (O. 0, O. 0, 333) 

WRITE  (*,412> 

READ   (*, 200)  ANSWER 

IF   ((ANSWER  .£Q.      '  Y'  )  .OR.   (ANSWER  . EQ.  '  y'  )  >  GOTO  15 

13   WRI""E(*,  417) 

READ (*, 200)   ANSWER 

IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' ) )  THEN 

****  FILL  0'  s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 
DO  106  I  =  1,31 
DO  106  J  =  1,31 
ZF(I,J)  =0.0 
106  CONTINUE 

ZFMAX  =  -3. 3E20 
ZFMIN  =  3. 3E20 
DN  =  (N-l ) /2. 0 
P  =  6. 283185 
DO  107  I  =  1,N 
DO  107  J  =  1,  N 
RLPART   =  0. 0 
IMGPART  =  0.0 
DO  108  L  =  1,N 
DO  108  K  =  1 , N 

Rl(l)  =  COS (-P* (L-l ) * ( I-DN-1 ) /N) 
Rl (2)  =  SIN <-P*(L-l )*( I-DN-1 ) /N) 
R2!l)  =  COS (-P*(K-1 ) *(J-DN-1 ) /N> 
R£(2)  =  SIN(-P*(K-1) *(J-DN-1) /N) 
RLPART   =  RLPART   +•  Z (L, K ) * ( Rl ( 1 ) *R2 ( 1 ) 

*  -Rl (2) *R2 (2) ) 
IMGPART  =  IMGPART  +  Z (L, K ) * ( Rl ( 1 ) *R2 (2) 

*  +R1 <£) *R2 ( 1 ) ) 
108        CONTINUE 

ZF(I,J)  =  SCRT(RLPART**2  +  IMGPART**2) 
IF  (ZF(I,J)   . GT.   ZF^AX;  'HEN 
ZFMAX  =  ZF(I, J) 
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107 


ENDIF 

IF  <ZF(I,J)  . LT.  ZFMIN)  THEN 

ZFMIN  =  ZF(I,J> 
ENDIF 
CONTINUE 


*****  OUTPUT  THE  ZF  MATRIX  ***** 
WRITE  (*, £05)  ' ***   FOURIER  TRANSFORMATION 
WRITE  (*,£1£) 
DO  103  I  =  1,N 

WRITE  (*,300>   (ZF(I,J),  J  =     1,N) 
WRITE  <*, £10) 
109  CONTINUE 

WRITE  (*,213) 


',N,  '  X  ',IM,  ' 


30 


WRITE  (*,418) 

READ   <*, £00)  ANSWER 

IF   ( (ANSWER  . NE.  ' Y' ) 


AND.   (ANSWER  .  NE.   '  y'  )  )  *30T0  16 


***** 

AS^ 

<  THE  PARAMETi 

WRITE 

(*, 

£10) 

WRITE 

(*. 

*)  '  ***   E  N 

WRITE 

(*, 

405) 

READ 

(* 

*)   AZIM 

WRITE 

(* 

406) 

READ 

(* 

*)   ELEV 

WRITE 

(* 

^08) 

READ 

(* 

*)    ITRIM 

WRITE 

(* 

403) 

READ 

(* 

, *)   IDIV 

WRITE 

(* 

411) 

READ 

(* 

, 133)   CTEXT 

WRITE 

(* 

,  40  1 ) 

READ 

(* 

,£00)   ANSWER 

E  T  E  R  S 


*♦*' 


*****  INITIALIZE  PL0T38  ***** 


OR.   (ANSWER  . EQ.  ' y' ) )   THEN 


IF   (  (ANSWER  . EQ.  ' Y'  ) 
CALL  PLOTS (0, 0, £) 

ELSE 

CALL  PLOTS ( 0, 33, 33 ) 

ENDIF 

WRITE  C*, 420) 

READ   (*, 200)   ANSWER 


CALL  WINDOW(XLOL,  YLOL,  XUPR, YUPR) 

IF   ((ANSWER  .  EQ.  ' Y'  )  .OR.   (ANSWER  . EQ.  ' y'  )  )   THEN 
DLEV  =  (ZFMAX-ZFMIN) /FLOAT (N) 
CALL  ZLEVEL (ZF, 31,31, N, N, DLEV, ZLEV, N+l  ) 
DO  110  I  =  1,N+1 
LD I G ( I  )  =  £ 
LWGT(I)  =  1 
110      CONTINUE 


£33  C 
£34 

£36 


*****  DRAW  THE  CONTOUR  MAP  ***** 

CALL  ZCNTUR(ZF, 31 , 3 1 , N, N, 0. 5, 0. 5, 7. 5, 5. 5, ZLEV, LDIG, LWGT, 

N-t-1,  0.  10,  10) 
CALL  SYMBOL <5. 5, 0. O, O. £,' CONTOUR  MAP', 0.0, 11) 
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Li  ne# 
£37 
£38 
239 
£40 
£41 
£4£ 
£43 
£44 
£45 
£46 
£47 
£48 
£49 
£50 
£51 
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ELSE 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 
CALL  MESHS(ZF,  31 ,  31 ,  N,  N,  AZ IM, ELEV, 0. 5, 0. 5,  7.  5, 5. 5, 

3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 
*****  ANNOTATION  OF  THE  GRAPH  ***** 
CALL  SYMBOL (5. 5, 0. 3, 0. £,' AZIMUTH:     ',0.0,10) 
CALL  NUMBER (999. 0, 999.0, 0. £, AZIM, 0. 0, £) 
CALL  SYMBOL (5. 5, 0. 0, 0. £, ' ELEVATION:' , O. 0, 10) 
CALL  NUMBER (999. O, 999. 0, 0. £,  ELEV, O. 0, £) 
DY  =  (ZF(1, 1) /90. O)  *  ELEV 
CALL  P3D£D(1.  0,  1.  0,  ZF(1,  D-DY,  XR,  YR) 
CALL  SYMBOL  (XR,  YR,  0.  £5,  '  *'  ,  0.  0,  1) 
CALL  SYMBOL ( 1. 0, 0. 1, 0. £,' *  =  ORIGIN' ,0. O, 10) 

ENDIF 

CALL  SYMBOL ( 1. O, S. 75, 0. £5, CTEXT, 0. 0, £0) 
CALL  SYMBOL  (6.  0,  6.  5,0.  £,  '  £-D  DFT',0.0,7) 
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£87 

£88 

289 
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291 


(ANSWER  .EQ.  ' y' ))  GOTO  30 


16 


*****  OUTPUT  THE  GRAPH  ***■< 

CALL  PLOT (0. 0, O. O, 999) 

WRITE  (*, 412) 

READ   (*, 200)  ANSWER 

IF   ((ANSWER  .EQ.  » Y*  )  .OR. 
ENDIF 

WRITE  (*,412) 
READ   (*,20O)  ANSWER 

IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' ) )  GOTO  10 
STOP 


199  FORMAT (A20) 

£00  FORMAT (A) 

£05  FORMAT (/, £0X, A£5, I£, A3, I£, A8, /) 

210  FORMAT (/) 

211  FORMAT (/,  A8,  12, A47) 

212  FORMAT(/, 2X, ' (AZIMUTH  320. 0) ' , 46X, ' ( AZ IMUTH  230.0)',/) 

213  FORMATf/, 2X, ' (AZIMUTH  050. 0) ', 46X, ' (AZ IMUTH  140.0)',/) 
300  FORMAT ( 10 (F7.  2,  IX)  ) 

399  FORMAT (/, 5X,  A4,  \) 

400  FORMAT (9X,  \) 

401  FORMAT (/, 5X, ' SEND  GRAPH  TO  THE  PRINTER(Y  or  N):  ',\) 

402  FORMAT (/, 5X, ' NUMBER  OF  ROWS/COLUMNS  FOR  R  AND  S(l  to  25):  ',\) 

403  F0RMAT(5X, A4,  12, A3,  \) 

404  FORMAT (5X,  A2,  12, A5,  \) 

405  FORMAT(/, 5X, ' AZIMUTH(0. 0  to  360.0  DEGREES):  '  ,\) 

406  FORMAT (/, 5X, ' ELEVATION (90. 0  to  -90.0  DEGREES):  ',\> 

408  FORMAT (/, 5X, ' TRIM (0=N0, l  =  Xs, 2=Ys) :  ',\) 

409  FORMAT (/, 5X, ' 2, 4  OR  8  SUBGRIDS:  ',\) 

411  FORMAT(/, 5X, ' TITLE  OF  GRAPH(UP  TO  20  CHAR):  ',\> 

412  FORMAT  (/,  5X,  '  DO  YOU  WANT  TO  CHANGE  PARAMETERS'7  ',\) 

413  FORMAT  (/,  5X,  '  DO  YOU  WANT  TO  REPEAT  THE  PROCESS'1  ',\> 

417  FORMAT(/, 5X, ' DO  YOU  WANT  FOURIER  TRANSFORMATION  ?  ',\) 

418  FORMAT (/, 5X, ' DO  YOU  WANT  TO  MAKE  GRAPH  ?  '  ,\) 

419  FORMAT (/, 5X, 'DO  YOU  WANT  TO  CORRECT  ?  ' , \ ) 

420  FORMAT (/, 5X, ' DO  YOU  WANT  CONTOUR  MAP  ?  ',\) 
END 
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Al 


Type 

REAL 


Offset  P  Class 
26 


127 


D  Line*  1 

A£ 

REAL 

30 

03 

REAL 

34 

A4 

REAL 

38 

ANSWER 

CHAR*1 

74 

PZIM 

REAL 

114 

Bl 

REAL 

4£ 

B2 

REAL 

46 

CI 

REAL 

50 

C£ 

REAL 

54 

COS 

INTRINSIC 

CTEXT 

CHRR*£0 

126 

DLEV 

REAL 

£16 

DN 

REAL 

166 

DY 

REAL 

146 

ELEV 

REAL 

na 

FLOAT 

INTRINSIC 

I 

INTEGER* 

£ 

60 

IDIV 

INTEGER* 

£ 

124 

IMGPAR 

REAL 

190 

IPROJ 

INTEGER* 

p 

£2 

ITRIM 

INTEGER* 

2 

122 

J 

INTEGER* 

36 

K 

INTEGER* 

£ 

202 

L 

INTEGER* 

194 

LDIG 

INTEGER* 

60O0 

LARGE 

LWGT 

INTEGER* 

6062 

LARGE 

MASK 

INTEGER* 

0 

LARGE 

N 

INTEGER* 

58 

NRNG 

INTEGER* 

£4 

P 

REAL 

170 

R 

REAL 

0 

LARGE 

Rl 

REAL 

0 

LARGE 

R£ 

REAL 

a 

LARGE 

RLPRRT 

REAL 

186 

S 

REAL 

£500 

LARGE 

SIN 

INTRINSIC 

SQRT 

INTRINSIC 

U 

REAL 

76 

VERTEX 

REAL 

0 

LARGE 

XLQL 

REAL 

£ 

XR 

REAL 

150 

XUPR 

REAL 

10 

YLC'L 

REAL 

6 

YR 

REAL 

154 

YUPR 

REAL 

14 

Z 

REAL 

5000 

LARGE 

ZF 

REAL 

8844 

LARGE 

ZFMAX 

REAL 

158 

ZFMIN 

REAL 

162 

ZLEV 

REAL 

12688 

LARGE 

ZLOW 

REAL 
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Name 


Type 


Size 


CI  ass 


MAIN 

MESHS 

NUMBER 


PROGRAM 

SUBROUTINE 

SUBROUTINE 


128 


D  Line* 

P3D20 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

ZCNTUR 

ZLEVEL 
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Pace 

03-25-33 

12:24:27 
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SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 


Pass  Ore 


No  Errors  Detected 
291  Source  Lines 


0> 


129 


APPENDIX     B  Jaqs  1 

03-56-95 

£1 : 09:36 

D  Line**  1.  7                                    Microsoft  FQRTR.AN77  V3.  20  02/3- 

1  3S70RAGE:   3 

2  SPAGE3IZE:53 

3  c  *********+**.**+»**********************************+**+***^ 

4  C  *                                                                       * 

5  C  *   THE  PURPOSE  OF  'HIS  PROGRAM  IS   TO   COMPUTE   AND  GRAPH  THE     * 
S  C  *   FREQUENCY  RESPONSE  OF  A  £-0  DIGITAL  FILTER.                       * 

7  C  *                                                                           + 

8  C  *                        EVANGELOS  THEQFILOU                           * 
g  C  ^*****^*********************»**********************»************ 

10  C  PROGRAM  2D-DATA-FIELD 
11 

12  C  *****  VARIABLE  DECLARATIONS  ***** 

13  REAL  A(7,  7) , 3(7, 7) ,  Rl  (7, 7, 2) ,  R2(7, 7, 2) , 

14  *  RLPART, IMGPART, Z<51, 51) , 

15  *  VERTEX ( 16) , ZLEV (51) 

16  INTEGER  MASK (3000) , LDIG ( 51 ) , LWGT ( 51 ) 

17  CHARACTER*!  ANSWER 
IB  CHARACTER*20  CTEXT 
19 

20  DATA  XLOL/0. 0/, YLOL/0. 0/, XUPR/3. 5/, YUPR/7. O/, 

21  *  ZLOW/1. 0E25/, IPROJ/0/, NRNG/100/ 
22 

23  C  *************   MAIN    PROGRAM   ************* 
24 

25  10   WRITE  (*, 401) 

26  READ   (*, *)     IT 

27  IF  (IT  . GT.  25)  GOTO  10 
33  WRITE  (*,  4-02) 

39  READ   (*, *)  K 

30  K  =  K  +  1 

31  WRITE (*,*)  '  ENTER  VALUES  OF  COEFFICIENTS:' 

22  DO  100  I  =  0,K-1 

1  33  DO  100  J  =  0, K-l 

2  24  WRITE  (*,404)  '  3  ('  ,  I,  '  ,  '  , J,  '  )  :  ' 
2     25  READ  (*,*)  .9(1  +  1,  J-t-1) 

2     26    100  CONTINUE 
27 

23  DO  101  I  =  0, K-l 

1  29  DO  101  J  =  0,K-1 

2  40  WRITE  (*,  404)  '  AC  ,  I,  '  ,  '  ,  J,  '  )  :  ' 
2     41  READ  (*,*)  A(I  +  1,J+1> 

2     42    101  CONTINUE 
43 

4.4  WRITE  (*,  419) 

45  READ   (*,200)  ANSWER 

46  IF   ((ANSWER  .ED.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' ) )  GOTO  10 

47 

ua    C  ****  FILL  0' s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 

49  DO  107  I  =  1,51 

1  50  DO  107  J  =  1,  51 

2  51  Z (I, J)  =  0. 0 
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107  CONTINUE 


54 
55 
56 
57 
58 
59 
60 
SI 
62 
S3 
64 
S5 
66 
S7 
68 
S3 
70 
71 
72 
73 
74 
75 
76 
77 
73 
79 
80 

ai 

82 
32 
34 
85 
86 
37 
38 
39 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
100 
101 
102 


=  3.  9E20 

=  -3.  9E20 

3.  14153 

=  E*P  /  CIT-1) 

-o    -    STEP 


=  U 


ZMIN 
ZMflX 

3    = 

STEP 

Wl  = 

L 

DO  102  I  =  1,  IT 
Wl  =  Wl  +  STEP 
W£  =  -P  -  STEP 
DO  103  J  =  1, IT 
L  =  L  +  1 
W2  =  W2  ->-  STEP 


104 


105 


103 
102 


DO  104  m  = 
DO  104  N 
Rl (M+l 
Rl  (M+l 
R2(M+1 
R2(M+1 

CONTINUE 

RLNOM 

IMGNOM 

RLDEN 

IMGDEN 

DO  105 
DO  105 
RLNOM 


0,K-1 
=  O.K-1 
N+l 
N+l 

N+l 
N+l 


1  ) 
2) 
1) 
2) 


=  COS(-M 
=  5IN(-M 
=  COS(-N 
=  SIN(-N 


Ull) 

Wl  ) 
W2> 
W£) 


=  0. 

=  0. 

=  0. 

=  0. 

"<1  = 


o 

0 

0 

0 

0,K-1 

=  0,K-1 

=  RLNOMn 


B  (M+l, Nh 


1 ) + (Rl  (M+l, N+l,  1  )  *R2 (M+l, N+i,  1  ) 

-  Rl (M+l, N+l, 2) *R2 (M+l, N+l, 2) ) 
I MGNOM+B  (  M+  1 ,  N+  1  )  *  (  R 1  (  M+  1 ,  N+  1 ,  1 ) *R2  <  M+ 1 ,  N+ 1 ,  2 ) 

+  R2(M+1, N+l, 1 ) *R1 (M+l, N+l, 2) ) 
RLDEN   =  RLDEN+A (M+l, N+l ) *(R1 (M+l, N+l, 1 ) *R2 (M+l, N+l, 1 ) 

-  Rl (M+l, N+l, 2) *R2 (M+l, N+l, 2) ) 
MGDEN+A (M+l, N+l ) * ( Rl (M+l, N+l, 1 ) *R2(M+1, N+l, 2) 

+  R2 (M+l, N+l, 1) *R1 (M+l, N+l, 2) ) 


IMGNOM  = 


IMGDEN  = 


CONTINUE 

ELEMENT  =  SORT (RLN0M**2 
SQRT (RLDEN**2 
ZCI, J)  =  ELEMENT 


IMGNCM**2> 
IMGDEN**2) 


/ 


C. 


IF  (Z(I,J)  .  GT.  ZMflX)  THEN 

ZMPX  =  Z(I,J) 
END  IF 
IF  (Z(I,J)  . LT-  ZMIN)  THEN 

ZMIN  =  Z (I, J) 
END  IF 
CONTINUE 
CONTINUE 

**♦**  OUTPUT  THE  Z  MATRIX  *♦*** 
WRITE  (*,  205)  '***■*■*   Z   M  A  T  R 
WRITE  (*,212) 


',  IT,  '  X  ',  r 


*■»*♦*' 
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110 
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114 
115 
116 
117 

ua 

113 

ISO 
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123 
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7 
DO  106  I  =  1,  IT 

WRITE  (*,200)   (Z(I,J),  J  =  1,  IT) 
WRITE  (*, 210) 
106  CONTINUE 

WRITE  (*,213) 

WRITE  <*,418> 

READ   (*, 200)  ANSWER 

IF   ((ANSWER  . NE.  ' Y' )  .AND.   (ANSWER  . NE.  ' y' ) )  GOTO  15 


*****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 
WRITE  (*, 210) 


WRITE 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 


108 


* )  ' **** 

4 1 0 ) 

*)   AZIM 

411) 

*)   ELEV 

413) 

*)    ITRIM 

414) 

*)    IDIV 

415) 

199)  CTEXT 
451) 

200)  ANSWER 


E  N  T 


PLOT   PARAMETERS 


*****  INITIALIZE  PLQT8S  ***** 


OR.   (ANSWER  . EQ.  ' y' ) )   THEN 


IF   ( (ANSWER  . EQ.  ' Y' ) 
CALL  PLOTS (O, 0, 2) 

ELSE 

CALL  PLOTS (0, 99,  99) 

END  IF 

WRITE  (*, 420) 

READ  (*,200)  ANSWER 


CALL  WINDOW (XLOL, YLQL, XUPR, YUPR) 

IF   ((ANSWER  .EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' ) )   THEN 
DLEV  =  (ZMAX-ZMIN) /FLOAT (IT) 
CALL  ZLEVEL (Z, 51,51, IT, IT, DLEV, ZLEV, IT+1 ) 
DO  108  I  =  1, IT+1 
LDIG( I)  =2 
LWGT(I)  =  1 
CONTINUE 

*****  DRAW  THE  CONTOUR  MAO  ***** 

CALL  ZCNTUR (Z, 51, 51,  IT,  IT, 0.  5, 0. 5, 8. 25, 6. 5,  ZLEV,  LDIG, LWGT, 
*  IT+1, 0. 10, 10) 

CALL  SYMBOL (5. 5, 0. 0, 0. 2, ' CONTOUR  MAP', 0.0, 11) 
ELSE 
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160 

161 
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164 
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168 

169 
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171 

172 
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174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

134 

185 

186 

137 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

138 

199 

200 

301 

202 

203 

204 


***** 

COLL 
y 

***** 
COLL 
COLL 
COLL 
COLL 
DY  = 
COLL 
COLL 
COLL 
END  I F 
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Microsoft  FGR7RAN77  V3.  20 
DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

MESHS (Z,  51,  51,  IT,  IT,  AZIM,  ELEV,  0.  5,  0.  5,  3.  25,  3.  5,  IDIV,  0 
3,  IPROJ,  1,  ZLGUI,  3,  ITRIM,  MASK,  VERTEX) 
ANNOTATION  OF  THE  GRAPH  ***** 

SYMBOL (5.  5,  0.  3,  0. 2,  ' AZIMUTH:     '  , 0. 0,  10) 

NUMBER (999. 0, 999. O, 0.  2, AZIM, 0. 0, 2) 

SYMBOL (5. 5, 0. 0, 0. 2,  ' ELEVATION: '  , 0. 0,  10) 

NUMBER  (999.  0,  999.  0,  0.  2,  ELEV,  0.  0,  2) 

(Z ( 1,  1 ) /30. 0)  *  ELEV 

P3D2D( 1. O,  1.  0,  Z (1,  1) -DY,  XR, YR) 

SYMBOL  (XR,  YR,  0.  25,  '  *'  ,  O.  O,  1) 

SYMBOL ( 1. O, 0. 1, 0. 2, ' *  =  ORIGIN' , 0. 0, 10) 


-22      - 

-22-35 
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02/34 


CALL  SYMBOL ( 1.  0,  3.  75,  0.  25,  CTEXT,  O.  0, 20) 

COLL  SYMBOL (3. 0, 3. 5, 0. 2, ' 2-D  DATA  FIELD' , 0. 0, 14) 

*****  OUTPUT  THE  GRAPH  ***** 
CALL  PLOT (0. 0, 0. 0, 393) 

WRITE  (*,41S) 

READ   (*,200)  ANSWER 

IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' > )  GOTO  20 

WRITE  (*,4i7) 

READ   (*, 200)  ANSWER 

IF   ((ANSWER  .EQ.   ' Y'  )  .OR.   (ANSWER  .20.  ' y'  )  )  GOTO  10 

STOP 
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205 
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3QQ 

4O0 

451 
401 
402 
404 
410 
41  1 
413 
414 
415 
413 
417 
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413 
420 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORmAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


( A20 ) 

(0) 

(/, 20X, A25, 12, A3, 12, A3, /) 

(  ) 


5X, A60) 


3lO. 0 ) 
050. 0) 


<*6X, 
45X, 


(AZIMUTH 
(AZIMUTH 


2X, ' (AZIMU'H 
2X, ' (AZIMUTH 
}<F7. 2,  IX)  ) 
X,  \) 
5X,'SEND  GRAPH  TO  THE  PRINTER (Y  or 
5X, ' DIMENSION  OF  OUTPUT  MATRIX (1  t 
5X,' ORDER  OF  TRANSFER  FUNCTION (0  t 
X,  A2,  II,  A,  II,  A3,  \) 
5X,  '  AZIMUTH(0. 0  to  330.0  DEGREES): 
5X, ' ELEVATION (30. 0  to  -90.0  DEGREE 
5X, 'TRIM (0=N0, l=Xs, 2=Ys) :  '  ,  \) 
5X,'2, 4  OR  3  3UBGRIDS:  ' , \) 
5X, ' TITLE  OF  GRAPH (UP  TO  20  CHAR): 
DO  YOU  WANT  TO  CHANGE  PARAMETE 
TO  REPEAT  THE  PROC 
TO  MAKE  GRAPH  ?  '  , 
TO  CORRECT  ?  ' , \) 
CONTOUR  MAP  ?  ' , \) 


£30 
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.  0  l  '  ,  /  ) 
.  0 )  '  ,  /  ) 


N)  :  '  ,  \) 
i  25)  :  '  ,  \) 


4) 

'  ,  \ 

S)  : 


) 


5X , 
5X, 
5X , 
5X, 
5X, 


DO 
DO 
DO 
DO 


YOU 
YOU 
YOU 
YOU 


WANT 
WANT 
WANT 
WONT 


RS  ? 

ESS 

\) 


) 

?  '  ,  \) 
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Name 


ype 


Offset  P  Clas'j 


fl 

REPL 

£ 

PNSWER 

CHAR*1 

i  a  1 1  o 

PZIM 

REAL 

18202 

B 

REPL 

138 

COS 

INTRINSIC 

CTEXT 

CHPR*20 

18214 

DLEV 

REAL 

18234 

DY 

REPL 

18244 

ELEMEN 

REPL 

18  130 

ELEV 

REOL 

13206 

FLOPT 

INTRINSIC 

I 

INTE5ER*£ 

18082 

IDIV 

INTEGER*£ 

13212 

IMGDEN 

INTEGER*2 

18176 

IMGNOM 

INTEGER*£ 

18170 

IMGPPR 

REAL 

****■»■ 

IPROJ 

INTEGER*£ 

18074 

IT 

INTEGER*£ 

18078 

ITRIM 

INTEGER*£ 

18210 

J 

INTEGER*£ 

18030 

K 

INTEGER*£ 

18080 

L 

INTEGER*£ 

18132 

LDIS 

INTEGER*£ 

17850 

LWGT 

INTEGER*£ 

1735£ 

M 

INTEGER*£ 

18150 

MP3K 

INTEGER*£ 

11850 

N 

INTEGER*£ 

18158 

NRNG 

INTEGER*£ 

18076 

P 

REAL 

1  3 1  £0 

Rl 

REPL 

334 

R£ 

REPL 

786 

RLDEN 

REPL 

18172 

RLNOM 

REPL 

18166 

RLPPRT 

REPL 

■**-»■■** 

SIN 

INTRINSIC 

SORT 

INTRINSIC 

STEP 

REPL 

18124 

VERTEX 

REPL 

11786 

Wl 

REPL 

18123 

W£ 

REPL 

18140 

XLOL 

REPL 

18054 

XR 

REPL 

18248 

XUPR 

REPL 

13062 

YLOL 

REPL 

18053 

YR 

REPL 

13352 

YUPR 

REPL 

1 8066 

Z 

REPL 

1178 
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D  Line**  1 
ZLEV    REAL 

ZLCW  RESL 
ZMAX  REAL 
ZMIN    REAL 


1  1  532 
18070 

lans 

iaii£ 


Micros* 


09-££-35 

r0RTRSN77  VS. £0  Oi/BH 


Name 


Type 


Size 


Class 


MAIN 
MESHS 

NUMBER 

P3D£D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

ZCNTUR 

ZLEV EL 


PROGRAM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 


Pass  One 


No  Errors  Detected 
£05  Source  Lines 
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U  ) 

O 
CD 

iT" 


i'; ; 

il 
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APPENDIX    C 


Line*     I  7 

1  SSTORAGE:   2 
£  SPQ6E5I2E:58 


""  icrosor  z 


Page  .  1 

03-26-35 

£0 : 4-2 : 32 
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7 

a 

3 
10 

n 

12 

13 
14 
15 
16 
17 

ia 

13 
20 
21 
22 
22 
24 
25 

as 

27 
29 
29 
30 

31 


34 


*****************+♦******♦****♦*♦♦♦-*-»■*-»■♦*-«.-».-•.*-«■-*-#. -«.-».*  *•»■■»-*  -* -*-«••«.  ■**-».  -«.-».♦♦ 

*  «■ 

*  THE  PURPOSE  OF  THIS  PROGRAM  IS   TO   COMPUTE   «ND  GRAPH  THE  * 

*  EQUATIONS  OF  ROBERT  P. R0ES3ER  IM  THE  "DISCRETE  STATE-SPACE  + 

*  MODEL  FOR  LINEAR  IMAGE  PROCESSING".   IT  TRANSFORMS  ALSO  THE  * 

*  OUTPUT  MATRIX  Y  ACCORDING  TO  FOURIER  ANALYSIS.  * 

*  * 

*  EVANGELOS  THEOFILOU  * 

PROGRAM  2D-DATA-FIELD 

**-♦*+  VARIABLE  DECLARATIONS  ***** 

REAL         Rl  (26, 26) ,  R2(£6,  26) , SI (26, 26) , S2(£6, 25) , 
»  FR1 (2) , FR2(2) , TRM(4, 4>  ,  IV  <4)  ,  OV  (4)  , IMGPART 

CHARACTERS  ANSWER 

****  VARIABLE  DECLARATIONS  FOR  PL0T33  ***** 
CHARACTER*£0  CTEXT 

COMMON        /WORK  /Z (26,26) , ZF(26, 26) , ZLEV(25) ,LDIG(£6) , 
*  LWGT ( £6 ) , MASK ( 2000 ) , VERTEX ( 1 £ ) 


DATA 


XLOL/0. 0/, YLOL/0. 0/, XUPR/3. 5/, YUPR/7. 0/, 
ZLCW/1. 0E25/, IPROJ/0/, NRNG/ 100/ 


MAIN 


PROGRAM 


**♦♦  ■»♦*-*•♦■*•»•*♦ 


*  ASK  THE  REQUIRED  VALUES  FOR  THE  MODEL 
E  (*,  403) 


10  WRI 

READ   (*, *)  KK 
IF  ( (KK  . LT.  3) 


.OR.   (KK  .  G" 


25) )  GOTO  10 


25 

DO  100  I  =  1,KK  +  1 

L 

36 

DO  100  J  =  1, KK+ 

3 

37 

Rl  (I,  J)  =  0.  0 

3 

33 

R2(I,J)  =  0.0 

3 

39 

SI  (I, J)  =  0. 0 

3 

40 

S2CI,  J)  =  0.  0 

3 

41 

100 

CONTINUE 

*2 

43 

DO  101  I  =  1,4 

1 

44 

DO  101  J  =  1,4 

3 

A5 

TRM ( I, J)  =  0 

3 

46 

101 

CONTINUE 

47 

43 

DO  102  I  =  1,4 

1 

49 

IV (  I)  =0. 0 

1 

50 

OV(I)  =0. 0 

1 

51 

102 

CONTINUE 
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52 

53  WRITE  <*,211)  'ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  Rl<#. #) ' 

54  DO  103  I  =  1,KK 

1     55  WRITE  (*,  404)  ' Rl < 1 , ' , I , ' ) :  ' 

1      56  READ   (*,  *)  Rl  (  1,  I  ) 
1     57    103  CONTINUE 
58 

59  WRITE  (*, 211)  'ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R2(#. #)' 

60  DO  104  I  =  1,KK 

1      61  WRITE  (*,  404)  '  R2 ( 1 ,  '  ,  I ,  '  )  :  ' 

1     62  READ   (*, *)  R2( 1, I) 
1      63    104  CONTINUE 
64 

55  WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  SI  <**.**>  ' 
66  DO  105  I  =  1,KK 

1     67  WRITE  (*, 405)  ' S 1  ( '  ,  I ,  '  ,  1 )  :  ' 

1     68  READ   (*,  *>  SKI,  1) 
1      69    105  CONTINUE 
70 

71  WRITE  (*, 211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S2(#. #)  ' 

72  DO  106  I  =  1,KK 

1     73  WRITE  (*, 405)  ' S2 ( ' , I , ' , 1 ) :  ' 

1      74  READ   (*, *>  S2( I, 1 ) 
1     75    106  CONTINUE 
76 

77  WRITE  (*, 211)  'ENTER  VALUES  FOR  THE  OUTPUT  VECTOR (#.#)  ' 

78  0V(1)  =  1 

79  WRITE  <*, 409)  'bOl:  ' 
60  READ   (*,  *)  OV (3) 

31  WRITE  (*,409)  'aOl:  ' 

82  READ   (*,*)  0V(4> 
83 

84  WRITE  (*, 211)  'ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX (#.  #)    ' 

85  TRM(1,2)  =  1 

86  TRM(4, 1)  =  1 

87  WRITE  (*,409)  'alO:  ' 

88  READ   (*, *)  TRM< 1, 1) 

89  WRITE  (*,409)  ' a20 :  ' 

90  READ   (*, *)  TRM<2, 1) 

91  WRITE  (*, 409)  ' bll :  ' 

92  READ   (*,  *)  TEMP 

33  TRM(1,3>  =  TEMP  +  OV  (  3 ) *TRM ( 1 ,  1 ) 

34  WRITE  (*, 409)  'all:  ' 

35  READ   (*,  *)  TEMP 

36  TRM(1,4)  =  TEMP  +  OV (4) *TRM < 1 , 1 ) 

37  WRITE  (*, 403)  'b21:  ' 

38  READ   (*,*)  TEMP 

33  TRM(£,3)  =  TEMP  +  OV ( 3 ) *TRM (2, 1 ) 

100  WRITE  (*,403)  'a21:  ' 

101  READ   (*,*)  TEMP 

102  TRM(2,4)  =  TEMP  +  OV (4) *TRM (2, 1 ) 
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111 
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1 
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2 

120 

2 
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2 
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2 

123 

a 
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2 
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126 
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123 

129 

130 
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133 
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1 
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1 
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1 
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139  C 

140 

1 

141 

D 
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£ 
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1 
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3 
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£ 
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1 
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1 
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20  :  <*2  :  Z: 
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TRM(4, 3)  =  0V (3) 
TRM (4, 4)  =  OV (4) 

WRITE  (»,  211)  'ENTER  VALUES  FOR  THE  INPUT  VECTOR  ( *.  I*)  ' 

IV<3>  =  1 

WRITE  <*, ^09)  ' bOO:  ' 

READ  <*,  *>  IV(4) 

WRITE  (*,409)  'blO:  ' 

READ  (*,  *)  TEMP 

IV<1>  =  TEMP  +  IV(4> *TRM(1,  1) 

IV(2)  =  IV (4) *TRM(2, 1) 

U  =  1.0 

DO  107  I  =  1,KK 

DO  107  J  =  1,KK 

R1(I+1,J)  =  TRM( 1, 1) *R1 ( I, J)  +  R2(I,J)  +  TRM ( 1 , 3) *S  1  (  I ,  J)  + 

*  TRM( 1, 4)*S2(I, J)  +  IV(1)*U 
R2(I  +  1,J>  =  TRM(2,  1)  *R1  (  I,  J)  +  TRM  (2.  3)  +S1  (  I ,  J*  ■+• 

*  TRM (2, 4) *S2( I, J)  +  IV(2)*U 
SI (I, J+l)  =  U 

S2(I,J+1)  =  Rl(I,J)  +  0V(3)*S1(I,J)  +•  OV  (4)  *S2(  I,  J)     +     IV(4)*U 
U  =  0.  0 

107  CONTINUE 

WRITE  (*,  205)  '**♦■*■*  INPUT  VECTOR  +****' 
WRITE  (*, 300)   (IV(I),I=  1,4) 

WRI~E  (*,205)  '*****  OUTPUT  VECTOR  *****' 
WRITE  C*, 300)   (OV(I),I  =  1,4) 

WRITE  (*, 205)   '  **■*•**  TRANSITION  MATRIX  *****' 
DO  108  I  =  1,4 

WRITE  (*,300)   (TRM(I,J),J  =  1,4) 

WRITE  (*,210> 

108  CONTINUE 

****  FILL  0' s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  *•**■* 
DO  109  I  =  1,26 
DO  109  J  =  1, 26 
HI,  J)     =0.0 

109  CONTINUE 

DO  110  I  =  1,KK 

DO  110  J  =  1, KK 

Z(I,J)  =  R1(I,J)  +  0V(3)*S1(I,J)  +  OV  (4)  *S2  tl,  J) 

110  CONTINUE 

WRITE     (*, 205)     '*****    Rl        MATRIX        '  ,  KK,  '      X     '  ,  KK,  '     ** *-*-><•' 
DO     111     I     =     1 , KK 

WRITE     (*, 300)      (Rid,  J),     J    =     1,KK) 

WRITE     (*,210) 
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111     CONTINUE 


WRITE     (*, 205)     '*****    R2       M    A    T    R     I     X       » , KK, '     X     '  ,  KK,  '     *****> 
DO     112     I     =     1,KK 

WRITE     <*, 300)      (R2(I,J>,     J    =     t,  KK) 

WRITE     (#,210) 

112  CONTINUE 

WRITE  (*,  205)  '*****  SI   M  A  T  R  I  X   '  ,  KK,  '  X  '  ,  KK,  '  *****' 
DO  113  I  =  1,KK 

WRITE  <♦, 300)   (SI (I, J),  J  =  1,KK) 

WRITE  (*,210) 

113  CONTINUE 

WRITE  (*,205)  '♦**♦#  S2   M  ft  T  R  I  X   ' , KK, '  X  ' , KK, '  *****' 
DO  114  I  =  1,KK 

WRITE  (*, 300)   (S2(I,J),  J  =  1,  KK) 

WRITE  (*, 210) 

114  CONTINUE 

*****  OUTPUT  THE  Y  MATRIX  ***** 

WRITE  (*, 205)  '*****   Z   M  ft  T  R  I  X   ' , KK, '  X  ' , KK, '   *****> 

WRITE  (*,212) 

DO  115  1=  1, KK 

WRITE  (#,300)   <Z(I,J),  J  =  1,KK) 

WRITE  (#,210) 

115  CONTINUE 
WRITE  (#,213) 

WRITE(*, 419) 

REftD   (*,  200)  ANSWER 

IF   ((ANSWER  . NE.  'Y')  . ftND.   (ANSWER  . NE.  ' y' ) )  GOTO  21 

*****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 

20   WRITE  (#,210) 

WRITE  (*,  *)  ' ****   ENTER   PLOT   PARAMETERS   ****' 

WRITE  (*, 410) 

READ  (*, *)   AZIM 

WRITE  (*, 41 1 ) 

READ  (*,*)   ELEV 

WRITE  (*,413) 

READ  (*, *)    ITRIM 

WRITE  (*, 414) 

REftD  (*,  *)    ID  IV 

WRITE  (*,415) 

REftD  (*, 199)   CTEXT 

WRITE  (*, 451) 

REftD  (*, 200)   ANSWER 


*****  INITIALIZE  PL0T88  ***** 

IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ. 


y'  )  )   THEN 
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£05 
£06 
207 
£08 
£09 
£10 
211 
212 
£13 
£14 
£15 
216 
£17 

£ia 

£  1 9 
££0 
££1 
ooo 

op  2 

£24 
££5 
££S 
££7 
££3 
£29 
£30 
£31 
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1       7 

CALL  PLOTS (0,0,  2) 

ELSE 

COLL  PLOTS (0, 99, 99) 
END  IF 

COLL  WINDOW ( XLOL, YLOL, XUPR, YUPR) 


*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 
COLL  MESHS(Z,£6,  £6,  KK,  KK,  AZIM,  ELEV,  0.  5,  0.5,  3.  £5,  6.5,  IDIV,  0, 
►  3, IPROJ, 1, ZLQW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (1. 0, 6. 75, 0. £5, CTEXT, 0. 0, £0) 

CALL  SYMBOL (6. 0, 6. 5, 0. £, ' £-D  DATA  FIELD' , 0. O, 14) 

CALL  SYMBOL (5. 5, 0. 3, 0. £,' AZIMUTH:     ',0.0,10) 

CALL  NUMBER (999.  0,  999. 0, 0. £, AZIM, 0. 0, £) 

CALL  SYMBOL (5. 5, 0. 0, O. £, ' ELEVATION: ' , 0. 0, 10) 

CALL  NUMBER (999. O, 999. O, 0. £, ELEV, O. O, £) 

DY  =  (Z  (1,  1) /90. 0)  *  ELEV 

CALL  P3D£D  (1.0,  1.  0,  Z  (  1,  1 )  -DY,  XR,  YR) 

CALL  SYMBOL (XR, YR, 0. £5,  '  *'  , 0.  0,  1) 

CALL  SYMBOL ( 1. O, 0. 1, O. £,' *  =  ORIGIN' , 0. 0, 10) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (0. 0, 0. 0, 999) 

WRITE  (*,416) 

READ   (*, 200)  ANSWER 

IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  .ED.   ' y' ) )  GOTO  £0 


£34 
£35 
£36 
£37 
£33 
£39 
£40 
£41 
£4£ 
£43 
£44 
£45 
£46 
247 
£48 
£49 
£50 
251 


£54 
£55 


£1    WRITEt*, 418) 

READ(*,£00)   ANSWER 
IF   (  (ANSWER  . EQ.  ' Y'  ) 


OR. 


(ANSWER  . EQ.  ' y' ) )  THEN 


****  FILL  0' S  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS 
DO  116  I  =  1,£6 
DO  116  J  =  1,26 
ZF(I,J)  =0.0 
116    CONTINUE 

ZFMAX  =  -9. 9E£0 
ZFMIN  =  9. 9E£0 
DK  =  (KK  -  1)  /  £. 0 
P  =  3. 141592 
DO  117  M  =  1, KK 
DO  117  N  =  1, KK 
RLPART   =  0. 0 
IMGPART  =  O. O 
DO  118  L  =  1,KK 
DO  118  K  =  1,  KK 

FR1(1)  =  COS (-£*P*(L-1 ) *(M-DK-1 ) /KK) 
FR1(£)  =  SIN(-£*P*(L-1) *(M-DK-1) /KK) 
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Lineft 

4 

25b 

4 

257 

4 

253 

4 

253 

4 

260 

4 

£61 

4 

262 

Ua 

263 

2 

£64 

2 

265 

£ 

£66 

£ 

267 

2 

268 

£ 

263 

•3 

£70 

£71 

£7£ 

£73 

£74 

£75 

1 

£76 

1 

277 

1 

£78 

£73 

£80 

£81 

£S£ 

283 

284 

285 

286 

287 

288 

£83 

230 

231 

232 

233 

£34 

235 

236 

£37 

238 

233 

300 

301 

302 

303 

304 

305 

306 

118 


117 


113 


30 


Microsoft  F0RTRAN7' 

FR£(1)  =  C0S(-2*P*(K-1 > *(N-DK-1 ) /KK) 
FR2(2)  =  SIN <-2*P* <K-1 )* (N-DK-1 ) /KK) 
RLPRRT   =  RLPART   +  Z ( L, K) * ( FR1 ( 1 ) *FR£ ( 1 ) 

-FR1 (2) *FR£(2> ) 
IMGPART  =  IMGPART  +  Z <L, K) * ( FR 1 < 1 ) *FR2 (2) 

+FR1 (2) *FR2< 1 ) ) 
CONTINUE 

ZF<M,N)  =  SQRT(RLPART**2  +  IMGPART**2) 
IF  <ZF(M,N)  . GT.  ZFMAX)  THEN 
ZFMAX  =  ZF(M,N) 
ENDIF 
IF  (ZF(M,N)  . LT.  ZFMIN)  THEN 

ZFMIN  =  ZF(M,N> 
ENDIF 
CONTINUE 


V3. 


Pac  = 

5 

03-26- 

-85 

20  :  42 

:32 

20  02 

'84 

*****  OUTPUT  THE  ZF  MATRIX  ***** 

WRITE  (*, 205)  ' ***  FOURIER  TRANSFORMATION 

WRITE  (*,212) 

DO  113  I  =  1,KK 

WRITE  (*,300)   (ZF(I,J),  J  =  1,KK) 

WRITE  <*, 210) 
CONTINUE 
WRITE  (*,213) 


',KK,'   X 


KK, 


WRITE(*, 413) 

READ   (*, 200)  ANSWER 

IF   (  (ANSWER  .  NE.  '  Y' 


***** 

ASK  THE  PARAMETERS  POR 

THE 

WRITE 

(*, 210) 

WRITE 

(*,  *)  '  **♦   ENTER 

P  L 

WRITE 

(*,  4  10) 

READ 

(*,  *)   AZIM 

WRITE 

(*, 411) 

READ 

(*,  *)   ELEV 

WRITE 

(♦,  4-13) 

READ 

(*,*)    ITRIM 

WRITE 

(*, 414) 

READ 

(*,*)    IDIV 

WRITE 

(*, 415) 

READ 

(*, 133)   CTEXT 

WRITE 

(», 451 ) 

READ 

(*,200)   ANSWER 

AND.   (ANSWER  . NE.   'y'))  GOTO  2£ 


GRAPH  ***** 


QT       PARAMETERS 


*****     INITIALIZE    PL0TS8    ***** 

IF        ((ANSWER    . EQ.      ' Y' )     .OR.      (ANSWER 

CALL  PLOTS (0,0,2) 
ELSE 

CALL  PLOTS (O,  33, 33) 
ENDIF 


EQ. 


)  ) 


THEN 
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307 

308  WRITE  C»,  420) 

309  READ  (*,£00)  ANSWER 
310 

311  COLL  WINDOW  (XLQL,  YLCL,  XLiPR,  YUPR) 

313  IF  ((ANSWER  . EQ.  ' Y')  .OR.   (ANSWER  .ED.  '  y'  ) )  THEN 

314  DLEV  =  (ZFMAX-ZFMIN) /FLOAT(KK) 

315  CALL  ZLEVELtZF, 26,26, KK, KK, DLEV, ZLEV, KK+1) 

316  DO  136  I  =  1,KK+1 
1    317  LDIG(I)  =  2 

1    318  LWGT(I)  =  1 

1    313  136      CONTINUE 

320  COLL  ZCNTUR(ZF,  26,  26,  KK,  KK,  0.  5,  0. 5, 3. 25, 6. 5,  ZLEV, LDIG, LWGT, 

321  *  KK+1, 0.  10,  10) 

322  COLL  SYMBOL  (5.  5,  O.  0,  0.  2,  '  CONTOUR  MOP',  0.0,  11) 

323  ELSE 

324-  C           *****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

325  CALL  MESHS(ZF,  26,  26,  KK,  KK,  AZIM, ELEV, 0. 5, 0. 5, 8. 25, 6. 5,  IDIV,  0, 

326  *  3,  IPROJ,  1, ZL0W,3,  ITRIM, MASK, VERTEX) 
327 

328  C          *****  ANNOTATION  OF  THE  GRAPH  ***** 

323  CALL  SYMBOL  (5.  5,  0.  3,  0.  2,  '  AZIMUTH:     ',0.0,10) 

330  CALL  NUMBER (333.  0,  933.  0,  0.  2, AZIM, O. 0, 2) 

331  CALL  SYMBOL (5. 5, 0. 0, 0. 2, ' ELEVATION:' , 0. 0, 10) 

332  CALL  NUMBER<333.  0,  333.  O,  0.  2,  ELEV, O. O, 2) 

333  DY  =  (ZF(1, 1) /30. 0)  *  ELEV 

334  CALL  P3D2DU.  0,  1.  0,  ZF(1,  1  ) -DY,  XR,  YR) 

335  CALL  SYMBOL  (  XR,  YR,  0.  25,  '*',  0.  0,  1) 

336  CALL  SYMBOL ( 1. 0, 0. 1, 0. 2, ' *  =  ORIGIN' , 0. 0, 10) 

337  ENDIF 

338  CALL  SYMBOL (1.  O, 6.  75,  0.25,  CTEXT, 0. 0, 20) 
333  CALL  SYMBOL  (6.  0,  S.  5,  0.  2,  '  2-D  DFT',0.0,7) 
340 

341  C  *****  OUTPUT  THE  GRAPH  ***** 

342  CALL  PLOT (0. 0, 0. 0, 333) 

343  WRITE  (*, 416) 

344  READ   (*, 200)  ANSWER 

345  IF   ((ANSWER  .ED.  ' Y' )  .OR.   (ANSWER  .ED.  ' y' ) )  GOTO  30 

346  22    ENDIF 

347  WRITE  (*, 417) 

348  READ   (*, 200)  ANSWER 

343  IF   ((ANSWER  .ED.  ' Y' )  .OR.   (ANSWER  .ED.  ' y' ) )  GOTO  10 

350  STOP 

351 

352  133  FORMAT (A20) 

352  200  FORMAT (A) 

354  205  FORMAT(/, 1SX, A23, 12, A3, 12, AS, /) 

255  210  FORMAT () 

356  211  FORMAT (/, 5X, A56) 

357  212  FORMAT (/, 2X, ' (AZIMUTH  320. 0) ' , 46X, ' ( AZ IMUTH  220.0)',/) 
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Line#  1 

7 

353 

213 

FORMAT (/, 2X, ' 

353 

300 

FORMAT ( 10 (F7. 

360 

400 

F0RMAT(9X, \> 

361 

403 

FORMAT (/, 5X, ' 

365 

404 

FORMAT (5X, A5, 

363 

405 

F0RMAT(5X, A3, 

364 

406 

FORMAT (5X, A2, 

365 

407 

FORMAT (5X, A2, 

366 

403 

FORMAT (5X, A3, 

367 

403 

FORMAT (5X, AS, 

363 

410 

FORMAT (/, 5X, ' 

363 

411 

FORMAT(/, 5X, ■ 

370 

413 

FORMAT (/, 5X, ' 

371 

414 

FORMAT (/, 5X, ' 

372 

415 

FORMAT (/, 5X, ' 

373 

416 

FORMAT</, 5X, ' 

374 

417 

FORMAT (/, 5X, ' 

375 

413 

FORMAT (/, 5X, ' 

376 

413 

FORMAT (/, 5X, ' 

377 

420 

FORMAT (/, 5X, ' 

373 

451 

FORMAT (/, 5X, ' 

373 

END 

Microsoft  FCRTRAN77  73 
(AZIMUTH  050. 0) ', 46X, ' (AZIMUTH  140.0)',/) 
2,  IX)  ) 

DIMENSION  OF  OUTPUT (K  = 1 to20 ) :  ' , \) 

12, A3, \) 

12, A5, \) 

12, A4, \) 

12,  12, A3,  \) 

12, A3, \) 

\) 

AZIMUTH (0.0  to  360.0  DEGREES):  ',\) 

ELEVATION (90. 0  to  -90.0  DEGREES):  ',\) 

TRIM(0=NO, l=Xs, 2=Ys) :  ',\) 

2, 4  OR  3  SUBGRIDS:  ' , \) 

TITLE  OF  GRAPH(UP  TO  20  CHAR):  ',\) 

DO  YOU  WANT  TO  CHANGE  PARAMETERS'1  '  ,\) 

DO  YOU  WANT  TO  REPEAT  THE  PROCESS"  '  ,\> 

DO  YOU  WANT  FOURIER  TRANSFORMATION  "?    ',\) 

DO  YOU  WANT  TO  MAKE  GRAPH  ?  ',\) 

DO  YOU  WANT  CONTOUR  MAP  ?  ',\) 

SEND  GRAPH  TO  THE  PRINTER(Y  or  N):  ',\) 
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Name 


Type 


Offset  P  Class 


ANSWER 

CHAR*1 

11060 

AZIM 

REAL 

11062 

COS 

INTRINSIC 

CTEXT 

CHAR*20 

11074 

DK 

REAL 

11114 

DLEV 

REAL 

11  163 

DY 

REAL 

1  1094 

ELEV 

REAL 

1 1 066 

FLOAT 

INTRINSIC 

FR1 

REAL 

10332 

FR2 

REAL 

1 0390 

I 

INTEGER*2 

1 0956 

IDIV 

INTEGER*2 

11072 

IMGPAR 

REAL 

11142 

IPROJ 

INTEGER*2 

10950 

ITRIM 

INTEGER*2 

11070 

IV 

REAL 

10393 

J 

INTEGER*2 

10964 

K 

INTEGER*2 

1  1  154 

KK 

INTEGER*2 

10954 

L 

INTEGER*2 

11146 

LDIG 

INTEGER*2 

5512 

/WORK   / 

LWGT 

INTEGER*2 

5554 

/WORK   / 

M 

INTEGER+2 

1  1122 

MASK 

INTEGER*2 

5616 

/WORK   / 

N 

INTEGER*2 

11130 

144 
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NRNG 

INTEGER* 

□V 

REAL 

P 

REPL 

Rl 

REPL 

R2 

REAL 

RLPPRT 

REPL 

SI 

REAL 

S2 

REPL 

SIN 

SQRT 

TEMP 

REAL 

TRM 

REAL 

U 

REAL 

VERTEX 

REPL 

XLOL 

REAL 

XR 

REPL 

XUPR 

REPL 

YLOL 

REAL 

YR 

REAL 

YUPR 

REAL 

Z 

REAL 

ZF 

REPL 

ZFMRX 

REAL 

ZFMIN 

REAL 

ZLEV 

REPL 

ZLQU 

REPL 

10952 

10314 

una 

3706 

11133 

5410 

8114 


10996 

i  oa  i  a 

1 1 000 
11616 
10930 

liosa 

10933 
10934 
11  102 
10942 
0 

£704 
11106 
11110 

5403 
10946 


INTRINSIC 
INTRINSIC 


/WORK 


/WORK 
/WORK 
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/WORK   / 


Name 


Type 


Size 


CI  ass 


MAIN 

MESHS 

NUMBER 

P3D2D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

WORK 

ZCNTUR 

ZLEVEL 


11630 


PROGRPM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

COMMON 

SUBROUTINE 

SUBROUTINE 


Pass  One 


No  Errors  Detected 
i79  Source  Lines 
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1  SSTORAGE:   2 

2  3PAGESIZE:58 


THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  COMPUTE  AND  GRPPH  THE 
EQUATIONS  OF  ROBERT  P. ROESSER  IN  THE  "DISCRETE  STPTE-3PPCE 
MODEL  FOR  LINEAR  IMAGE  PROCESSING".  IT  TRANSFORMS  ALSO  THE 
OUTPUT  MATRIX  Y  ACCORDING  TO  FOURIER  ANALYSIS. 


4 

C 

5 

c 

6 

c 

7 

2 

a 

c 

3 

c 

10 

c 

11 

c 

12 

c 

12 

c 

14 

15 

c 

16 

17 

18 

13 

20 

c 

21 

22 

22 

24 

23 

26 

27 

23 

c 

23 

20 

c 

21 

jC 

23 

24 

25 

.25 

27 

28 

23 

40 

41 

42 

42 

44 

45 

46 

47 

48 

43 

50 

51 

*  EVANGELOS  THEOFILOU 

PROGRAM  2D-DATA-FIELD 


*****  VARIABLE  DECLARATIONS  ***** 

REAL  R(2S,25, 4) , 31 (26,26, 4) , 32(26, 2S  .  4)  , 

*  Rl (2) , R2(2) , TRM( 12, 12), IV(12),0V<12), IMGPART 
CHARACTERS  ANSWER 

****  VARIABLE  DECLARATIONS  FOR  PL0T88  ***** 

CHARACTER*20  CTEXT 

COMMON        /WORK  /Z (26, 26) , ZF(26, 26) , ZLEV(26) , LDIG(26) , 

*  LWGT ( 26 ) , MASK ( 2000 )  ,  VERTEX (IS) 

DATA  XLOL/0. 0/, YLOL/0. 0/, XUPR/8. 5/, YUPR/7. 0/, 

*  ZLOW/1. 0E25/, IPROJ/0/, NRNG/100/ 


MP.IN     PROGRAM 


*****  psk  THE  REQUIRED  VALUES  FOR  THE  MODEL 
10  WRITE  (*, 401) 

READ   (*, *)  N 

IF  (  (N  .  LT.   1)  .OR.   (N  . GT.  4))  GOTO  10 
2   WRITE  (*, 402) 

READ   (*, *)  M 

IF  (  (M  .  LT.   1)  .OR.   (M  . GT.  4))  GOTO  2 
2   WRITE  (*, 402) 

READ   (*, *)  KK 

IF  (KK  .  GT.  25)  GOTO  2 

DO  100  I  =  1,KK+1 
DO  100  J  =  1,  KK-f  1 
DO  100  L  =  1,  N 
R(  I,  J,  L)  =  0. O 
SI  (  I,  J,  L)  =  O.  0 
S2(I,  J,  L)  =  0.  O 
100  CONTINUE 

DO  101  I  =  l,N-t-2*M 
DO  101  J  =  1,N+£*M 
TRMCI, J)  =  0. 0 
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2 

52 

53 

1 

54 

1 

55 

1 

56 

57 

58 

59 

1 

60 

2 

61 

2 

62 

2 

63 

64 

65 

66 

1 

67 

a 

68 

2 

63 

2 

70 

71 

72 

73 

1 

74 

2 

75 

2 

75 

2 

77 

78 

73 

80 

81 

1 

82 

1 

83 

1 

84 

1 

85 

86 

1 

87 

1 

88 

1 

83 

1 

30 

31 

1 

32 

1 

33 

34 

1 

35 

1 

36 

37 

38 

39 

1 

100 

1 

101 

1 

102 
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101  CONTINUE 

DO  102  I  =  1,N+2*M 
IV ( I)  =  0. 0 
OV ( I )  =0.  0 

102  CONTINUE 

WRITE  <*, 311)  'ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R (».#)' 
DO  103  I  =  1,KK 
DO  103  J  =  1, N 

WRITE  (*,  404)  '  R'  ,  J,  Ml,  '  ,  I,  '  )  :  ' 

READ   (*,  *)  R(l,  I, J) 

103  CONTINUE 

WRITE  <*, 311)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  SI (#. #)  ' 
DO  104  I  =  1,KK 
DO  104  J  =  1,M 

WRITE  (*,  405)  '  31  ('  ,  J,  '  )  ('  ,  I,  '  ,  1)  :  ' 
READ   (*,  *)  SI ( I,  1, J) 

104  CONTINUE 

WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S2<#.#)  ' 
DO  105  I  =  1,KK 
DO  105  J  =  1 , M 

WRITE  (*,405)  '  S2C  ,  J,  '  )  ('  ,  I,  '  ,  1)  :  ' 
READ   (», *)  S2( I, 1, J) 

105  CONTINUE 

WRITE  (♦,311)  'ENTER  VALUES  FOR  THE  INPUT  VECTOR(#.#>  ' 

IV  (1)  =  1.0 
DO  106  I  =  1 , M 

WRITE  (*, 408)  'a(0',I,'):  ' 

READ  (*,*)  IV(N+I) 

TRM(N+I,N+1)  =  -IV(N-t-I) 

106  CONTINUE 

DO  107  I  =  1,M 

WRITE  (*,408)  'b(0',I,'):  ' 
READ  (*,*)   IV(N+M+I) 
TRM<N+M+I,  N+l )  =  -IV(N+M+I) 

107  CONTINUE 

DO  108  I  =  1,M-1 

TRM (N+I, N+I+l )  =  1.0 

108  CONTINUE 

DO  103  I  =  1,M-1 

TRM (N+M+I, N+M+I+l )  =1.0 
103  CONTINUE 

WRITE  (*,211)  'ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX  (4*.  #>    ' 
DO  110  I  =  1 , N 

WRITE  (*,406)  'a(',I,'0>:  ' 

READ   (*, *)  TEMP 

TRM( 1,1)  =  -TEMP 
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1 

103 

104 

105 

1 

106 

1 

107 
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1 

109 

2 

110 

^ 

111 

2 

112 

2 

113 

114 

1 

115 

2 

116 

£ 

117 

2 

118 

2 

113 

120 

1£1 

12£ 

123 

124 

125 

126 

1 

127 

1 

128 

1 

123 

1 

130 

131 

13£ 

133 

1 

134 

2 

135 

3 

136 

3 

137 

3 

138 

4 

139 

4 

140 

4 

141 

4 

142 

4 

143 

4 

144 

4 

145 

3 

146 

3 

147 

3 

148 

3 

143 

3 

150 

4 

151 

4 

152 

4 

153 
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110  CONTINUE 

TRM  (  1,  N+l  )  =  -1.0 
DO  111  I  =  a, ^ 

TRM(  1,1-1)  =  1.0 

111  CONTINUE 

DO  112  I  =  1,M 
DO  112  J  =  1, N 

WRITE  <*, 407)  'a(',J,  I,'):  ' 

READ   <*,  *)  TEMPI 

TRM(I+N,J)  =  TEMPI  +    TRM(1,J)  *  IV(N-t-I) 

112  CONTINUE 

DO  113  I  =  1,M 
DO  113  J  =  1, N 

WRITE  (*,407)  'b(',J,  I,'):  ' 

READ   (*,*)   TEMPI 

TRM(  I+N+M,  J)  =  TEMPI  +  TRM(1,J)  *  I V  <  N-t-M+I  ) 

113  CONTINUE 

WRITE  (*,211)  'ENTER  VQLUES  FOR  THE  OUTPUT  VECTOR (#.#)         ' 

WRITE  (*, 403)  ' b(OO) :  ' 

READ   (*,*)  TEMP 

0V(N+1)  =  -TEMP 

0V<N+M+1)  =  1.0 

DO  114  I  =  1,N 

WRITE  (*, 406)  ' b(' , I, ' O) :  ' 

REP.D   (*,*)  TEMPI 

OV  (  I )  =  TEMPI  +  TRM(1,I)  *  TEMP 

114  CONTINUE 

U  =  1.0 

DO  115  I  =  1,KK 

DO  115  J  =  1,  KK 

DO  116  II  =  1,N+2*M 

IF   (II  .  LE.  N>   THEN 
DO  117  J J  =  1,N+2*M 

IF   (JJ  . LE.  N)   THEN 

RU  +  1,  J,  II)=R(I  +  1,  J,  II)+TRM(II, JJ)*R(I, J, JJ) 
END  IF 
IF   ((JJ  .GT.  N>  .AND.   (JJ  . LE.  N+M))   THEN 

R(I+1, J, II)=R(I+1, J, II)+TRM( II, JJ) *S1 (I, J, JJ-N) 
END  IF 
117  CONTINUE 

R(I+1, J, II)=R(I+1, J, II)  +  IV(II)  *  U 
ENDIF 

IF   ((II  .GT.  N)  .  OND.   (II  . LE.  N+M))   THEN 
DO  118  JJ  =  1,N+2*M 

IF   (JJ  . LE.  N)   THEN 

SI  (I,  J-t-1,  II-N)  =  SI  (I,  J+l,  II-N)  -t-  TRM(II,JJ)  * 
*  R(I,J,  JJ) 
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154 

4 

155 

4 

156 

4 

157 

4 

153 

4 

159 

3 

160 

3 

161 

3 

162 

3 

163 

3 

164 

4 

165 

4 

166 

4 

167 

4 

163 

4 

163 

4 

170 

4 

171 

4 

173 

4 

173 

4 

174 

4 

175 

4 

176 

4 

177 

3 

173 

3 

179 

3 

ISO 

3 

131 

d 

133 

133 

134 

135 

136 

137 

133 

139 

190 

191 

1 

193 

1 

193 

1 

194 

195 

196  C 

197 

1 

193 

o 

199 

3 

200 

301 

303 

1 

303 

2 

304 
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END  IF 
IF  (<JJ  .ST.  NJ.flND. (JJ  . LE.  N+M))  THEN 

SI ( I, J+l, II-N)  =  SI ( I, J+l, II-N)  +  TRM ( II, JJ) * 

*  SI ( I, J, JJ-N) 
END  IF 

113  CONTINUE 

SI  (I, J+l,  II-N)  =  51  (I, J+l,  II-N)  +  IV(II)  *  U 
END  IF 

IF  (II  . GT.  N+M)  THEN 
DO  113  J J  =  1, N+£*M 

IF   (JJ  .  LE.  N)   THEN 

S3 (I, J+l, II-N-M)  =  S3 ( I, J+l, II-N-M)  +  TRM(II.JJ) 

*  *  R<I, J, JJ) 
END  IF 

IF    (  (JJ  .  GT.   N)   .AND.   (JJ  . LE.   N+M))    THEN 

S2(  I,  J+l,  II-N-M)  =  S3( I,  J+l,  II-N-M)  +  TRM(II,JJ) 

*  »  "SI  (  I,  J,  JJ-N) 
ENDIF 

IF   (JJ  .GT.  N+M)   THEN 

S£(I, J+l, II-N-M)  =  S£(I, J+l, II-N-M)  +  TRM(II,JJ) 

*  *  S3 ( I, J, JJ-N-M) 
ENDIF 

113  CONTINUE 

S3 (I, J+l, II-N-M)  =  S3 (I, J+l, II-N-M)  +  IV(II)  *  U 
ENDIF 
116      CONTINUE 

U  =  0.  0 
115  CONTINUE 

URITE  (*,205)  '****♦  INPUT  VECTOR  *****' 
WRITE  (*,300)   (IV (I), I  =  1,N+3*M) 

WRITE  (*,305)  ******  OUTPUT  VECTOR  *****' 
WRITE  (*,300)   (OV(I),I  =  1,N+2*M) 

WRITE  (*,305)  '****♦  TRPNSITION  MPTRIX  *■****' 
DO  130  I  =  1,N+3*M 

WRITE  (*, 300)   (TRM (I, J), J  =  1,N+3*M) 

WRITE  (*, 210) 

130  CONTINUE 

**«.*  FILL  O'  s  THE  TWO  DIMENTIONPL  GRID  OF  CONTROL  POINTS  *+*+ 
DO  131  I  =  1,36 
DO  131  J  =  1, 36 
Z(I,J)  =0.  0 

131  CONTINUE 

DO     133    I    =     1,KK 
DO    133    J    =     1,KK 

DO    133    LL    =    1, N+3*M 
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3    305  IF  (LL  . L£.  N)  THEN 

3    £06  Z(I,J)  =  Z(I,J)  +  OV(LL)  *  R(I.J,LL) 

3    £07  END  IF 

3    £08  IF  ( (LL  . GT.  N) . AND. (LL  . LE.  N+M)  )  THEN 

3    £09  Z(I,J)  =  Z(I,J)  +■  OV(LL)  *  SI  (I,  J,  LL-N) 

3    £10  END  IF 

3    £11  IF  (LL  . GT.  N+M)   THEN 

3    £1£  Z(I,J)  =  Z(I,J)  +  OV(LL)  *  S£(I,  J,  LL-N-M) 

3    £13  ENDIF 

3    £14  1£3      CONTINUE 

£    £15  1££  CONTINUE 

216 

£17  WRITE     (*,£05)     '*****    R       MATRIX        '  ,  KK,  '      X     »  ,  KK,  '     *****' 

£18        DO  124  I  ■  1,KK 
1    £19  DO  123  L  =  1, N 

£    ££0  WRITE  <*,  300)   (R(I,J,L),  J  =  1,KK) 

£    ££1  125    CONTINUE 
1    £££  WRITE  (*,210) 

1    ££3  124  CONTINUE 

££4 

££5  WRITE     (*,  £05)     '*****    SI        MATRIX        '  ,  KK,  '      X     '  ,  KK,  '     *****> 

££S        DO  12S  I  =  1,KK 

1  ££7  DO  127  L  =  1,M 

2  ££8  WRITE  <*,  300)   (SKI,  J,  L),  J  =  1,KK) 
£    ££9  127    CONTINUE 

1    £30  WRITE  (*, £10) 

1    £31  126  CONTINUE 

£33  WRITE     (#,£05)     '*****    3£       MATRIX        '  ,  KK,  '      X     '  ,  KK,  '      *****' 

£34         DO  128  I  =  1,KK 
1    £35  DO  129  L  =  1, M 

£    £36  WRITE  (*,  300)   (S£(I,J,L),  J  =  1,KK) 

£    £37  129    CONTINUE 
1    £38  WRITE  (*,£10) 

1    £39  128  CONTINUE 

£40 

£41  C      **•*■»*  OUTPUT  THE  Z  MATRIX  ***** 

£4£         WRITE  (*, 205)  '*****   Z   MATRIX   '  ,  KK,  '  X  '  ,  KK,  '    *****' 

£43        WRITE  (*,212) 

£44         DO  130  I  =  1,KK 
1    £45  WRITE  (*, 300)   (Z(I,J),  J  =  1,KK) 

1    £46  WRITE  (*, 210) 

1    £47  130  CONTINUE 

£48         WRITE  (*,£13) 

£49 

£50         WRITE(*,  419) 

£51         READ   (*,200)  ANSWER 

£5£         IF   ((ANSWER  . NE.  ' Y»  )  .AND.   (ANSWER  . NE.  ' y'  )  >  GOTO  £1 

£54  C      *****  OSK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 

£55  £0   WRITE  (*, £10) 
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£56 
£57 
£58 
£53 
£60 
£61 
£6£ 
£63 
£64 
£65 


£67 

£53 

£69 

£70  C 

£71 

£7£ 

£73 

£74 

£75 

£76 

£77 

£78 

£73  C 

£80 

£81 

£3£ 

£83  C 

£84 

£85 

£86 

£37 

£88 

£83 

£30 

£31 

£3£ 

£33 

£94 

£35  C 

£36 

£37 

£38 

£33 

300 

301   , 

302 

303 

304 

305  C 

306 


7 

WRITE  (* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (> 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  <* 

WRITE  (* 

READ  (* 

WRITE  (* 

READ  (* 


,  *)  '  **** 
,  410) 

,  *)   azim 

,411) 

, *)   ELEV 

,  413) 

, *)    ITRIM 

,  414) 

,*)    IDIV 

,  415) 

, 199)   CTEXT 

,451) 

,£00)   ANSWER 
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**♦♦' 


*****  INITIALIZE  PL0T88  ***** 


IF   ( (ANSWER  . EQ.  ' Y' ) 
CALL  PLOTS tO, 0, £) 

ELSE 

CALL  PLOTS (0,99,  99) 

ENDIF 


CALL  WINDOW(XLOL, YLOL, XUPR, YUPR) 


OR.   (ANSWER  . EQ.  ' y' ) )   THEN 


*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 
CALL  MESHS(Z,  £6.  £6,  KK,  KK,  AZIM,  ELEV,  0.  5,  0.  5,  3.  £5,  6.  5,  IDIV,  0. 
►  3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (1. 0,6.  75, 0. £5, CTEXT, O. 0, £0) 

CALL  SYMBOL (6. 0,6. 5, O. £, ' £-D  DATA  FIELD' , 0. O, 14) 

CALL  SYMBOL (5. 5, 0. 3, 0. £,' AZIMUTH:     ',0.0,10) 

CALL  NUMBER (999. O, 999. 0, 0. £, AZIM, O. O, £) 

CALL  SYMBOL (5. 5. 0. 0, 0.  £,  ' ELEVATION:'  , 0. 0,  10) 

CALL.  NUMBER  (999.  O,  999.  0,  0.  £,  ELEV,  0.  0,  £) 

DY  =  (Z ( 1, 1 ) /90. 0)  *  ELEV 

CALL  P3D£D ( 1. 0,  1. O,  Z  (  1,  1 ) -DY,  XR,  YR) 

CALL  SYMBOL (XR, YR, 0. £5, ' *' , 0. 0, 1) 

CALL  SYMBOL ( 1. O, 0. 1, O. £,' *  =  ORIGIN' , 0. O, 10) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (0. 0, 0.0, 999) 

WRITE  (*,41£) 

READ   (*, £00)  ANSWER 

IF   ( (ANSWER  . EQ.  ' Y' ) 


WRITEt*, 418) 
READ(*,£00)   ANSWER 
IF   ( (ANSWER  . EQ.  ' Y' ) 


OR.   (ANSWER  .EQ.  ' y' ) )  GOTO  £0 


OR.   (ANSWER  .EQ.  ' y' ) )  THEN 


****  FILL  0'  s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **-** 
DO  132  I  =  1,26 
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2 

309 
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312 

313 

314 
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1 

316 

2 
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2 

318 

2 

319 

2 

320 

4 

321 

4 

322 

4 

323 

4 

324 

4 

325 

4 

326 

4 

327 

4 

323 

4 

329 

5 

330 

p 

331 

132 


134 


334 

336 

336 

339  C 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352  C 
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DO  132  J  =  1,26 
ZF(I,J)  =0.0 
CONTINUE 

ZFMAX  =  -9. 9E20 
ZFMIN  =  9. 9E20 
DK  =  (KK  -  1)  /  2.0 
P  =  3. 141592 
DO  133  M  =  1, KK 
DO  133  N  =  1,KK 
RLPART   =  0.  0 
IMGPART  =  0. 0 
DO  134  L  =  1, KK 
DO  134  K  =  1,KK 
Rl  (1)  = 
Rl (2)  = 
R2(l)  = 
R2(2)  = 
RLPART 


COS ( -2*P* ( L- 1 ) * ( M-DK- 
S I N ( -2*P* (L-l ) *  < M-DK- 
COS <-2*P* (K-l ) * (N-DK- 
SIN (-2*P* (K-l ) *(N-DK- 
=  RLPART   +•  Z  (L,  K)  *(l 


-1 
-1 
-1 
-1 
Rl 
-Rl 
IMGPART  =  IMGPQRT  +  Z(L,K)*(R1 

+  R1 


)  /KK 

)  /KK 
)  /KK 
)  /KK 
(  1  )* 
(2)  * 
(  1  ) 
(2)  * 


) 
) 
) 

) 

R2  (  1 ) 

R2 (2) ) 

R2(2) 

R2( 1)  ) 


CONTINUE 

ZF(M,  IM)  =  SQRT <RLPART*+2  +  IMGPART**2) 

IF  (ZF(M,i\l>  . GT.   ZFMAX)  THEN 

ZFMAX  ■  ZF(M,N) 
END  IF 

IF  (ZF(M,N)  . LT.  ZFMIN)  THEN 
.   ZFMIN  =  ZF (M, N) 
END  IF 
CONTINUE 

*****  OUTPUT  THE  ZF  MATRIX  ***** 

WRITE  (*, 205)  ' ***  FOURIER  TRANSFORMATION   ' , KK, '  X  ' , KK, '  ***' 

WRITE  (*,212) 

DO  135  I  =  1,KK 

WRITE  (*,300)   (ZF(I,J),  J  =  1,  KK) 

WRITE  (*, 210) 
CONTINUE 
WRITE  (*,213) 

WRITE(*, 419) 

READ        (*, 200)     ANSWER 

IF        ((ANSWER    . NE.      'Y')     .AND.      (ANSWER    . NE.      ' y' ) )     GOTO    22 

*****    ASK    THE    PARAMETERS    FOR    THE    GRAPH    ***** 
WRITE     (*,210) 


WRITE  (* 

WRITE  (* 

READ  (* 

WRITE  (* 


*>     -  ***       ENTER       PLOT       PARAMETERS 

410) 

*)        AZIM 

411) 
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258  READ   (*, *)   £L£V 

359  WRITE  <*,413) 

360  READ   (*, *)   I TRIM 

361  WRITE  (*, 414) 

36£  READ   (*, *)    IDIV 

363  WRITE  (*,  4.15) 

364  READ   (*,  199)   CTEXT 

365  WRITE  (*, 451) 

366  READ   (*, 200)   ANSWER 
367 

368  C  *****  INITIALIZE  PL0T88  ***** 

369  IF   ((ANSWER  . EQ.  ' Y'  )  .OR.   (ANSWER  . EQ.  ' y'  )  )   THEN 

370  CALL  PLOTS  (0,0,2) 

371  ELSE 

372  CALL  PLOTS (O, 99, 99) 

373  ENDIF 
374 

375  WRITE  (*,  420) 

376  READ  (*, 200)  ANSWER 
377 

378  CALL  WINDOW ( XLOL, YLOL, XUPR, YUPR) 
379 

380  IF  ((ANSWER  .  EQ.  '  Y»  )  .OR.   (ANSWER  .  EQ.  '  y'  )  )  THEN 

381  DLEV  =  (ZFMAX-ZFMIN) /FLOAT (KK) 

382  CALL  ZLEVEHZF,  £6,  26,  KK,  KK,  DLEV,  ZLEV,  KK+1) 

383  DO  136  I  =  l.KK+1 
1    384  LDIG(I)  =2 

1    385  LWGT(I)  =  1 

1    386    136  CONTINUE 

387  CALL  ZCNTUR(ZF, 26,  £6, KK, KK,  O. 5, O. 5,8. £5, 6. 5,  ZLEV, LDIG, LWGT, 

388  *  KK+1, 0. 10, 10) 

389  CALL  SYMBOL (5. 5, 0. 0, 0. 2, 'CONTOUR  MAP' , O. O, 1 1 ) 

390  ELSE 

391  C  *****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

392  CALL  MESHS(ZF,  £6,  26,  KK,  KK,  AZIM, ELEV, 0. 5, 0.5, 8.25, 6. 5,  IDIV, 0, 

393  *  3,  IPROJ,  1,  ZLOW, 3,  ITRIM, MASK, VERTEX) 
394 

395  C  *****  ANNOTATION  OF  THE  GRAPH  ***** 

396  CALL  SYMBOL  (5.  5,  0.  3,  0.  2,  '  AZIMUTH:     ',0.0,10) 

397  CALL  NUMBER  (999.0,  999.  0,  0.2,  AZIM,  0.  0,  2) 

398  CALL  SYMBOL (5. 5, 0. 0, 0. 2, ' ELEVATION:' , 0. 0, 10) 

399  CALL  NUMBER (999. 0, 999. 0, 0. 2, ELEV, 0. 0, 2) 

400  DY  =  (ZF(1, 1) /90. 0)  *  ELEV 

401  CALL  P3D£D( 1. O, 1. 0, ZF (1, 1 ) -DY, XR, YR) 

402  CALL  SYMBOL ( XR, YR, 0. 25, '*', 0. 0, 1) 

403  CALL  SYMBOL ( 1. 0, 0. 1, O. 2, ' *  =  ORIGIN' , 0. 0, 10) 

404  ENDIF 

405  CALL  SYMBOL ( 1. 0, 6. 75, 0. 25, CTEXT, O. O, £0) 

406  CALL  SYMBOL  (6.  0,  6.  5,  0.  2,  '  2-D  DFT',  0.0,7) 
407 

408  C  *****  OUTPUT  THE  GRAPH  ***** 
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CALL  PLOT  (0.  0,  0.  0,  333) 

WRITE  (*,  416) 

READ   (*,2G0)  ANSWER 

IF   ((ANSWER  . EQ.  ' V )  .OR.   (ANSWER  . EQ.  ' y' ) )  GOTO  30 
END  IF 
WRITE  <*,417) 
READ   (*, 200)  ANSWER 

IF   ((ANSWER  . EQ.  ' Y'  )  .OR.   (ANSWER  .ED.  ' y'  )  )  GOTO  10 
STOP 


FORMAT (A£0) 

FORMAT (A) 

FORMAT (/, 18X, A£3, 12, A3, I2,fl8, /) 

FORMAT ( ) 

FORMAT(/, 5X, A46) 

FORMAT (/, £X, ' (AZIMUTH  3£0. 0) ' , 46X, 

(AZIMUTH  050.0)',46X, 

£,  IX)  ) 


FORMAT(/, £X, 
FORMAT ( 10 (F7 
FORMAT (3X, \) 
FORMAT(/, 5X, 
FORMAT(/, 5X, 
FORMAT (/, 5X, 
FORMAT (/, 5X, 
F0RMAT(5X, Al 
FORMAT (5X, A3 
FORMAT (5X, A£ 
FORMAT  <5X, A£ 
FORMAT (5X, A3 
FORMAT (5X, AS 
FORMAT(/, 5X, 
FORMAT(/, 5X, 
FORMAT?/, 5X, 
FORMAT(/, 5X, 
FORMAT(/, 5X, 
FORMAT (/, 5X, 
FORMAT(/, 5X, 
FORMAT (/, 5X, 
FORMATt/, 5X, 
FORMAT(/, 5X, 
END 


(AZIMUTH 
(AZIMUTH 


£30. 0) 
140. 0) 


/) 
/  ) 


\) 


SEND  GRAPH  TO  THE  3RINTER(Y  or  N):  ',\) 

NUMBER  OF  HORIZONTAL  STATES (N=l to4) 

NUMBER  OF  VERTICAL  STATES (M=lto4) :  ',\) 

DIMENSION  OF  OUTPUT ( lto£5) :  '  ,\) 

12, A3, 12, A3, \) 

I£, A£,  I£,  A5,  \> 

I£,  A4,  \) 

12,  :£, A3,\) 

12, A3, \) 

\  ) 

AZIMUTH(0. O  to  360.0  DEGREES):  ', 

ELEVATION (30. 0  to  -30.0  DEGREES): 

TRIM(0=NO, l=Xs, £=Ys) :  ',\) 

£, 4  OR  3  SU3GRIDS:  ' , \) 

TITLE  OF  GRAPH (UP  TO  £0  CHAR):  ', 

DO  YOU  WANT  TO  CHANGE  PARAMETERS'1 

DO  YOU  WANT  TO  REPEAT  THE  PROCESS 

DO  YOU  WANT  FOURIER  TRANSFORMATION  ?  '  ,\) 

DO  YOU  WANT  TO  MAKE  GRAPH  ?  »,\) 

DO  YOU  WANT  CONTOUR  MAP  ?  '  ,\) 


\) 


\) 


\) 


\) 
,  \) 


Name 


Type 


Offset  P  Class 


ANSWER  CHAR*1 

AZIM  REAL 
COS 

CTEXT  CHAR*£0 

DK  REAL 

DLEV  REAL   . 

DY  REAL 

ELEV  REAL 


33434 
33436 

33448 
33488 
33536 
33468 
33440 


INTRINSIC 
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Fi_3AT 

I 

INTEGER*2 

IDIV 

INTEGER+2 

II 

INTEGER*2 

IMGPPR 

REAL 

IPROJ 

INTEGER*2 

ITRIM 

INTEGER+2 

IV 

REAL 

J 

INTEGER*£ 

JJ 

INTEGER*2 

K 

INTEGER*^ 

KK 

INTEGER*2 

L 

INTEGER*S 

LDIG 

INTEGER*^ 

LL 

INTEGER*a 

LWGT 

INTEGER*^ 

M 

INTEGER*^ 

MASK 

INTEGER*^ 

IM 

INTEGER*2 

NRNG 

INTEGER*2 

OV 

REAL 

P 

REAL 

R 

REAL 

Rl 

REAL 

R2 

REAL 

RLPP.RT 

REAL 

SI 

REAL 

S2 

REAL 

SIN 

SORT 

TEMP 

REAL 

TEMPI 

REAL 

TRM 

REAL 

U 

REAL 

VERTEX 

REAL 

XLOL 

REAL 

XR 

REAL 

XUPR 

REAL 

YLDL 

REAL 

YR 

REAL 

YUPR 

REOL 

Z 

REAL 

ZF 

REAL 

ZFMOX 

REAL 

ZFMIN 

REAL 

ZLEV 

REAL 

ZLQW 

REAL 

33iea 

33446 
33336 
33512 
33158 
33444. 
33042 
33176 
33344 

33166 
33134 

5512 
33334 

556-+ 
33164 

5616 
33162 
33160 
33090 
33432 

33026 
33034 
33503 
10813 
£1634. 


33276 
33238 
32450 
33320 
11516 
33138 
33472 
33146 
33142 
33476 
33150 
0 

2704 
33480 
33434 

5408 
33154 
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INTRINSIC 


/WORK   / 


/WORK 


/WORK   / 


INTRINSIC 
INTRINSIC 


/WORK   / 


/WORK 
/WORK 


/WORK   / 


155 


D  Line* 


Page   11 

09-86-aS 

19 :3c:  :0b 

Microsoft  FGRTRAN77  V3.  £0  02/84 


Name 


Type 


Size 


Clas 


MO  IN 

ME5HS 

NUMBER 

P3D£D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

WORK 

ZCNTUR 

ZLEVEL 


usao 


PROGRAM 

SUBROUT 

SUBROUT 

SUBROUT 

SUBROUT 

SUBROUT 

SUBROUT 

SUBROUT 

COMMON 

SUBROUT 

SUBROUT 


INE 
INE 

INE 
INE 
INE 
INE 
INE 

INE 

INE 


Pass  One 


No  Errors  Detected 
448  Source  Lines 
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1  SLARGE 

2  SSTORAGE:   2 

3  *PAGESIZE:58 
4 

C      *  * 

*  THE  PURPOSE  OF  THIS  PROGRAM  13   TO   CODE  THE  1-D  (DISCRETE     * 

*  TIME)  SYSTEM  TO  A  £-D  SPACIAL  SYSTEM.  * 


5  C 

S  C 

7  C 

3  C 

9  C 

10  C 

11  c 

12  C 

IS 

14  C 

IS 

IS 

17 

18 

19  C 

SO 

El 

22 

23 

24 

25 

26 

27  C 

23 

23  C 

30 

31 

32 

33 

34 

35 

36 

37 

33 

33 

40 

41  C 

42 

1 

43 

z> 

44 

O 

45 

46 

47  C 

^8 

1 

43 

1 

SO 

1 

51 

*  EVANGELOS  THEOFILOU 

PROGRAM  2D-DATA-PIELD 


*****  VARIABLE  DECLARATIONS  ***** 

REAL         R  (25,  625)  ,  S  (25,  525)  ,  Rl  (2)  ,  R2  (2)  ,  TRM  (50,  50)  ,  IV  (50)  , 
h  IMGPART 

CHARACTER* 1  ANSWER 

***♦  VARIABLE  DECLARATIONS  FOR  PLCT33  ***** 

CHARACTER*20  CTEXT 

COMMON        /WORK  /Z (26, 26) , ZF(25, 26) , X (620) , V(630) , ZL£V(26) , 

*  LDIG(26) , LWGT(26) , MASK (3000) , VERTEX (IS) 

DATA  XLOL/0. 0/, YLOL/0. 0/, XUPR/3. 5/, YUPR/7. 0/, 

*  ZLOW/1. 0E35/, IPROJ/0/, NRNG/100/ 


♦♦♦♦**♦***♦*♦   MAIN    PROGRAM   *■**■#■**•*■*■*■»•*■*■* 
*****  flSK  THE  REQUIRED  VALUES  FOR  THE  MODEL  ***** 

2  WRITE  (*,  403) 
READ   (*, *)  N 

IF  (  (N  .LT.  3)  .OR.   (N  .GT.  23))  GOTO  2 

3  WRITE  (*,  *04.) 
READ   (*, *)  M 

IF  (  (M  . LT.  2)  .OR.   (M  .  GT.  25))  GOTO  3 

WRITE  (*, 401 ) 

READ  (*,200)  ANSWER 

IF   ((ANSWER  . EQ.  '  Y1  )  .OR.   (ANSWER  .  EQ.  '  y'  )  )  THEN 

♦***  FILL  0'  5  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  *♦** 
DO  36  I  =  1, 25 
DO  35  J  =  1, 25 
Z(I,J)  =0.0 
36    CONTINUE 


ENTER  VALUES  FOR  Y  MATRIX  **-< 
DO  37  I  =  1 , M*N 

WRITE  (*,  403)  '  YC  ,  I,  '  )  :  ' 
READ  C»,  *)  R(  1,  I) 
37    CONTINUE 
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52  END  IF 

53  IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  'y'))   GOTO  4 
54 

55  WRITE  (*, 402) 

56  READ  (*, 200)  ANSWER 

57  IF   ((ANSWER  .  EQ.  'V')  .OR.   (ANSWER  . EQ.  '  y'  >  )  THEN 

58  C  ****  INITIALIZE  THE  TRANSITION  MATRIX  **** 

59  DO  98  I  =  1,  M+N 

1  £0  DO  98  J  =  1 ,  M+N 

2  61  TRM(I,J)  =0.0 
2     62     98    CONTINUE 

63 

64  C  ***♦  ENTER  VALUES  FOR  TRANSITION  MATRIX  **** 

65  DO  99  I  =  1,M+N 

1  66  DO  99  J  =  1,  M+N 

2  67  WRITE(*,407)  ' T ( '  ,  I ,  '  ,  '  , J,  '  )  :  ' 
2     68  READ  (*,  *>  TRM(I,J) 

2     69     99    CONTINUE 

70  END  IF 
71 

72  C  *****  INITIALIZE  R  AND  S  ARRAYS  ***** 

73  DO  100  I  =  1,25 

1  74  DO  100  J  =  1, 625 

2  75  R(I, J)  =  0. O 
2  76  S(  I,  J)  =  0.  0 
2     77    100  CONTINUE 

73 

73  C  *****  INITIALIZE  INPUT  VECTOR  ***** 

80  DO  101  I  =  1,50 

1      81  IV(I)  -  0.0 
1      82    101  CONTINUE 

83 

84  WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R#' 

85  DO  102  I  =  1,M 

1     86  WRITE  (*,405)  '  R',I,':  ' 

1      87  READ   (*,  *)  R  (  I,  1) 
1      88    102  CONTINUE 
83 

90  WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S#' 

91  DO  103  I  =  1,N 

1      92  WRITE  (*,  405)  'S',I,':  ' 

1      93  READ   (*,  *)  S (  I,  1  ) 
1      94    103  CONTINUE 
95 

96  WRITE  (*,211)  'ENTER  VALUES  FOR  THE  INPUT  VECTOR' 

97  IV  (1)  =  1.0 

98  WRITE  (*, 406)  'aOl:  ' 

99  READ  (*,  *)  IV(M-t-l) 
100 

101  IF   ((ANSWER  .EQ.  ' Y'  )  .OR.   (ANSWER  . EQ.  '  y'  )  )   GOTO  5 

102  C  *****  INITIALIZE  TRANSITION  MATRIX  ***** 
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103  DO  104  I  =  1,35 

1  104  DO  104  J  =  1,25 

2  105  TRM ( I, J)  =  0.  0 
2    106  104  CONTINUE 

107  WRITE  <*,211)  'ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX' 

108  TRM  (M-t-1,  M-t-N)  =  -IV(M-t-l) 

109  WRITE  <*, 406)  'alO:  ' 

110  READ   (*,  *)  TEMP 

111  TRM(1,M)  =  -TEMP 

112  TRM(1,M+N)  =  -1.0 

113  WRITE  (*,406)  'all:  ' 

114  REP.D  (*,  *)  TEMP 

115  TRM(M+1,M)  =  TEMP  -  TRM(1,M)  *  TRM ( M+ 1 , M+N ) 
116 

117  DO  105  I  =  2, M 

1     1  13  TRM(  1,1-1)  =  1.0 

1     119  105  CONTINUE 

120 

121  DO  106  I  =  2+M,  M-t-N 

1     122  TRM  (I,  1-1)  =1.0 

1    123  106  CONTINUE 

124 

125  5   U  =  1. 0 

126  DO  107  I  =  1,N*M 

1  127  DO  108  J  =  1, M+N 

2  128  IF  (J  .  LE.  M)  THEN 

2  129  DO  109  JJ  =  l,M+N 

3  130  IF  (JJ  .  LE.  M)  R(J,  I-t-1)  =  R(J, 1+1)  + 
3    131  *                                 R ( JJ, I) *TRM (J; JJ) 
3    132  IF  (JJ  .  GT.  M)  R  (J,  I-t-1)  =  R(J,  I-t-1)  + 
3    133  *                              S (JJ-M, I) *TRM( J, JJ) 
3    134  109        CONTINUE 

2    135  R(J, 1+1)  =  R<J, 1+1)  +  IV<J)*U 

2    136  END  IF 

2    137 

2    138  IF  (J  .GT.  M)  THEN 

2  139  DO  110  JJ  =  1,M+N 

3  140  IF  (JJ  .LE.  M)  5(J-M,  I-t-1)  =  S(J-M,  1  +  1)  + 
3    141  *                                     R (JJ,  I) *TRM (J, JJ) 
3    142  IF  (JJ  .  GT.  M)  S(J-M,  I-t-1)  =  S(J-M,  I-t-1)  + 
3    143  *                                   S (JJ-M, I ) *TRM ( J, JJ) 
3    144  110         CONTINUE 

2     145  S(J-M,  I  +- 1  )  =  S(J-M,  1  +  1)  -t-  IV(J)*U 

2     146  ENDIF 

2    147  108    CONTINUE 

1     148  U  =  0. O 

1     149  107  CONTINUE 
150 

151  WRITE  (*,211)  '***-»•*  INPUT  VECTOR  ***-*•*' 

152  WRITE  (*,300)   (IV(I),I  =  l,M-t-N) 
153 
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WRITE  (*,  211)  '*****  TRANSITION  MATRIX  *****' 
DO  111  I  =  1 ,  M+N 

WRITE  (*,  300)   (TRM(I,J),J  =  1,M+N) 

WRITE  (*, £10) 

111  CONTINUE 

WRITE     (*, 311)     '*****    HORIZONTAL       STATES       R    *****> 
DO     112    I    =    1,M*N 

WRITE     (*, 300)      (R(J, I),     J    =     1,M) 

112  CONTINUE 

WRITE     (*, 211)     '        *****    VERTICAL      STATES       S    *****' 
DO    113    I    =    1,M*N 

WRITE     (*, 300)      (S(J, I),     J    =    1,N) 

113  CONTINUE 

****  FILL  0' s  THE  TWO  DIMENTI ONAL  GRID  OF  CONTROL  POINTS  **** 
DO  114  I  =  1,36 
DO  114  J  =  1,  £6 
2  (I,  J)  =0.0 

114  CONTINUE 

4   DO  US  I  =  1,11 

DO  115  J  =  1,  N 

Z.d,  J)  =  R(1,(I-1)  *N+J) 

115  CONTINUE 

*****  OUTPUT  THE  Y  ARRAY  ***** 
DO  113  I  =  1,M*N 

WRITE  (*, *>  R( 1, I) 
113  CONTINUE 

*****  OUTPUT  THE  Z  MATRIX  ***** 

WRITE     (*, £05)     '*****       Z       MATRIX        '  , M,  '      X     '  , N,  '        *****' 

WRITE  <*,£!£) 

DO  116  I  =  1,M 

WRITE  (*,300)   (Z(I,J),  J  =  1,N) 

WRITE  (*,£10) 

116  CONTINUE 
WRITE  (*,£13) 

WRITE(*, 421) 

READ   (*,£00)  ANSWER 

IF   ((ANSWER  .NE.  ' Y' )  .AND.   (ANSWER  . N£.  ' y' ) )  GOTO  13 

DO  117  I  =  1,630 

X(I)  =0.0 

Y ( I )  =0. 0 

117  CONTINUE 

DO  113  I  =  1,M*N 
X  (  I  )  =  I  *  1.0 
Y(I)  =  R(  1,  I) 
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CONTINUE 

WRITE  (*, 

415) 

READ   (*, 

139) 

CTEXT 

WRITE  (*, 

431) 

READ   (*, 

200 ) 

ANSWER 
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1    £05    113 

£06 

£07    18 

£08 

£09 

£10 

£11 

212  C  *****  INITIALIZE  PL0T88  ***** 

213  IF   ((ANSWER  . EQ.  ' Y' )  .OR.   (ANSWER  . EQ.  ' y' ) )   THEN 

214  CALL  PL0TS(O,O,2) 

215  ELSE 

£15  CALL  PLOTS (0,99, 99) 

£17  ENDIF 
£18 

£19  CALL  PLOT  (1.0,  1.0, -3) 

££0  CALL  SCALE (X, S. 0, M*N, 1 ) 

221  CALL  SCALE <Y, 4. O, M*N, 1) 

£22  CALL  STAXIS(0. 20, O. £0, 0.  Ill, 0.  112,  1) 

£23  CALL  AXIS (O. 0, 0. O,  '  X  AX  IS'  ,  -6,  6. 0, O. O, X (M*!\H-1 ) , X (M*N+2) ) 

224  CALL  AXISfO. 0, 0. 0,  ' Y  AX  IS'  ,  6,  4. 0, 90. O, Y (M*N+1 ) , Y (M*N+£)  ) 

££5  CALL  LINE(X, Y, M*N, 1, 0, 0) 

££6  CALL  PLOT(0. 0, 0. 0, -3) 

£27  CALL  SYMBOL ( 1. O, S.  75, 0. 25, CTEXT, 0. 0, 20) 

223  CALL  SYMBOL (6. 0, S. 5, 0. £,  '  1-D  DATA  FIELD'  , 0. O,  14) 
££9 

£30  C  *****  OUTPUT  THE  GRAPH  ***** 

£31  CALL  PLOT (0. 0, 0. 0, 399) 

£32  WRITE  (*,416) 

£33  READ   (*,  £00)  ANSWER 

£34  IF   ((ANSWER  .  EQ.  '  Y'  )  .OR.   (ANSWER  . EQ.  ' y'  )  )  GOTO  18 
£35 

£36    13   WRITE(*,419) 

£37  READ   (*, 200)  ANSWER 

£33  IF   ((ANSWER  .  NE.  '  Y'  )  .AND.   (ANSWER  . NE.  ' y'  )  )  GOTO  £1 
£33 

£40  C  *****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 
£41    £0   WRITE  (*,£10) 

242  WRITE  (*,*)  ' ****   ENTER   PLOT   PARAMETERS 

£43  WRITE  (*, 410) 

£44  READ   <*,*)   AZIM 

£4-5  WRITE  (♦,411) 

£46  READ   (*,  *)   ELEV 

£47  WRITE  (*, 413) 

£48  READ   (*, *)    ITRIM 

£43  WRITE  (*, 414) 

£50  READ   (*,  *)    ID  IV 

£51  WRITE  (*,415) 

£5£  READ   (*, 139)   CTEXT 

£53  WRITE  (*,451) 

£54  READ   (*, 20O)   ANSWER 
255 
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*****  INITIALIZE  PL0T38  ***** 

IF   ((ANSWER  .ED.  'V')  .OR.   (ANSWER  . EO.  ' y' ) )   THEN 

COLL  PLOTS (0,  0,  £) 
ELSE 

CALL  PLOTS (0, 33,  33) 
ENDIF 

CPLL  WINDOW (XLOL, YLCL, XUPR, YUPR) 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 
CALL  MESHS(Z,  £6,  £6,  N, M, AZIM, ELEV,  0.  5,  0.  5,  3. £5, 6. 5,  IDIV, 0, 
*  3,  IPROJ,  1,  ZLOW, 3,  ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (1.  0,  6.  75,  0.  £5, CTEXT, 0. 0, £0) 

CALL  SYMBOL (6. 0, 6. 5, 0. £,  ' £-D  DATA  FIELD'  , 0. O,  14) 


' AZIMUTH:     ' , 0. 0, 10) 
0.  £,  AZIM,  0.  O,  £> 
' ELEVATION:' , 0. 0, 10) 
O.  £,  ELEV,  0.  0,  £) 


CALL  SYMBOL (5. 5, O. 3, 0. £ 

CALL  NUMBER (333.  0,  333.  0 

CALL  SYMBOL (5. 5, 0. 0, 0.  £ 

CALL  NUMBER (333. 0, 333. 0 

DY  =  (Z  (1,  1) /30.  0)  *  ELEV 

CALL  P3D£D ( 1 . O,  1 .  0,  Z ( 1 ,  1 ) -DY,  XR, YR) 

CALL  SYMBOL ( XR, YR, 0. £5, ' *' , O. 0, 1 ) 

CALL  SYMBOL( 1. O, 0.  1, 0. £,  '  *  =  OR  I  GIN'  , O. O,  10) 


*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (0. 0, O. 0, 333) 

WRITE  (*,41S) 

READ   (*, £00)  ANSWER 

IF   ((ANSWER  . EQ.   '  Y'  )  .OR.   (ANSWER  . EO. 


y' ) )  GOTO  £0 


£1    WRITE(*,413) 

READ(*,£00)   ANSWER 

IF   ((ANSWER  .ED.  ' Y' )  .OR.   (ANSWER  .ED.  ' y' ) )  THEN 
;         ****  FILL  0' s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 
DO  13£  I  =  1,26 
DO  132  J  =  1,26 
ZF(I,J)  =0.0 
13£    CONTINUE 

ZFMAX  =  -3. 3E£0 
ZFMIN  =  3. 3E£0 
DN  =  (N-l ) /£. 0 
DM  =  (M-l) /£. O 
P  =  3. 14153E 
DO  133  MM  =  1,M 
DO  133  NN  =  1, N 

RLPART   =  0. O 

IMGPART  =  0. 0 

DO  134  L  =  1, M 
DO  134  K  =  1,N 
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Rl(l)  =  CaS(-£*P*<i_-l>*<MM-DI*-l)  /M> 
Rl (2)  =  SIN(-2*P*(L-1)*(MM-DM-1) /M) 
R2(l)  =  COS <-2*P* (K-l ) * (NN-DN-1 ) /N) 
R2(2)  =  SIN (-2*P* (K-l ) * (NN-DN-1 ) /N) 
RLPART   =  RLPART   -r-  Z (L, K) *(R1 ( 1 ) *R2 ( 1 ) 

-Rl (2) *R2(2) ) 
IMGPART  =  IMGPART  +  Z  (L,  !<>*(R1  (1  )*R2(2) 

+R1 (2) *R2 ( 1 ) ) 
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CONTINUE 
ZF(MM,NN)  = 
IF  (ZF(MM,NN> 
IF  (ZF(MM,NN) 
CONTINUE 


SORT (RLPART**2 
GT.  ZFMAX) 
LT.  ZFMIN) 


+  IMGPART**2) 
ZFMAX  =  ZF(MM,NN) 
ZFMIN  =  ZF (MM, NN) 


*****  OUTPUT  THE  ZF 

WRITE  (*, 205)  ' *** 

WRITE  (*,212) 

DO  135  I  =  1,M 

WRITE  (*,300>   (Z 
WRITE  <*,210) 

CONTINUE 

WRITE  (*,213) 

WRITE(*, 413) 


MATRIX 

FOURIER 


***** 
TRPNSFORMATION 


(I,  J)  ,   J  =  1,  N) 


READ 

(*, 200)  ANSWER 

IF   ( 

ANSWER  . NE.   '  Y'  )   .AND. 

(ANSWER  . NE-  ' y' ) 

)  GOTO  22 

***** 

ASK  THE  PARAMETERS  FOR 

THE  GRAPH  ***** 

WRITE 

(*, 210) 

WRITE 

(*,  *)  '  **-*   ENTER 

PLOT   p  p  R  p 

METERS   *** 

WRITE 

(*, 4  10) 

READ 

(*, *)   AZIM 

WRITE 

(*, 411) 

READ 

(*,  *>   ELEV 

WRITE 

(*, 413) 

READ 

(*,*)    ITRIM 

WRITE 

(*,  414) 

READ 

(*,*)    IDIV 

WRITE 

(*, 415) 

READ 

(*, 199)   CTEXT 

WRITE 

(*, 451 ) 

REPD 

<*,200)   ANSWER 

• 

*****  INITIALIZE  PL0T38  ***** 

IF   ((ANSWER  . EQ.  ' Y* )  .OR.   (ANSWER 

CALL  PLOTS (0, 0, 2) 
ELSE 

CALL  PLOTS (O,  99, 99) 
END  IF 

WRITE  (*,420) 


)  ) 


THEN 
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353  READ  (*, £00)  ANSWER 

353 

360  CALL  WINDOW (XLOL, YLOL,  XUPR, VUPR) 
361 

363  IF  ((ANSWER  .ED.   ' Y' )  .OR.   (ANSWER  .ED.  '  y'  )  )  THEN 

363  DLEV  =  (ZFMAX-ZFMIN) /FLOAT(M) 

364  CALL  ZLEVEL(ZF, 36,36, M,  N,  DLEV,  ZLEV, N) 

365  DO  136  I  =  1,N 
1    366  LDIG(I)  =  3 

1    367  LWGT(I)  =  1 

1    368  136      CONTINUE 

369  CALL  ZCNTUR(ZF,  36,  35,  M,  N,  0.  5,  0.  5,  8.  35,  6.  5,  ZLEV,  i_DIG,  LWGT, 

370  *  N,  0.  10,  10) 

371  CALL  SYMBOL  (5.  5,  0.  0,  O.  5,  '  CONTOUR  MAP'  ,  0.  0,  1  1  ) 
37£  ELSE 

373  C  *****    DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

374  CALl_  MESHS  (  ZF,  36,  36,  M,  N,  AZIM,  ELEV,  0.  5,  0.  5,  3.  35,  6.  5,  IDIV,  0, 

375  *  3,  IPROJ,  1,  ZLOW, 3,  ITRIM,  MASK,  VERTEX) 
376 

377  C  *****  ANNOTATION  OF  THE  GRAPH  ***** 

378  CALL  SYMBOL (5. 5, 0. 3, O. 3, ' AZIMUTH:    ',0.0,10) 
373  CALL  NUMBER (393.  0,  399.  0,  0.  3,  AZIM, 0. 0, 3) 

380  CALL  SYMBOL (5. 5, 0. 0, 0. 3, ' ELEVATION: ', 0. O, 10) 

381  CALL  NUMBER (999.  O,  999.  O,  0.  3, ELEV, 0. 0, 3) 
383  DY  =  (ZF(1, 1) /90. 0)  *  ELEV 

383  CALL  P3D3D(  1.  0,  1.  0,  ZF  (  1,  1  ) -DY,  XR,  YR) 

384  CALL  SYMBOL (XR,  YR,  0. 35,  '  *'  ,  0. 0,  1) 

385  CALL  SYMBOL  (  1.  0,  0.  1,  O.  3,  '  *  =  ORIGIN'  ,  O.  O,  10) 

386  ENDIF 

387  CALL  SYMBOL ( 1. 0, S.  75,  0.  35,  CTEXT, 0. 0, 30) 

388  CALL  SYMBOL  (6.  0,  6.  5,  0.  3,  '  3-D  DFT',0.0,7) 
389 

330  C  *****  OUTPUT  THE  GRAPH  ***** 

331  CALL  PLOT (0. 0, 0. 0, 333) 
333  WRITE  (*, 416) 

333  READ   (*,300)  ANSWER 

334  IF   ((ANSWER  .ED.   ' Y' >  .OR.   (ANSWER  .ED.  'y'))  GOTO  30 

335  33     ENDIF 

336  WRITE  (*, 417) 

337  READ   (■*,  200)  ANSWER 

338  IF   ((ANSWER  .ED.  ' Y' )  .OR.   (ANSWER  .ED.   ' y'  > )  GOTO  3 
333  STOP 

400 

401  133  FORMAT (A30) 

402  200  FORMAT (A) 

403  205  FORMAT (/,  18X, A23,  12, A3,  12, A8, /) 

404  210  FORMAT () 

405  211  FORMAT (/, 5X,  SOA) 

406  212  FORMAT (/, 2X, ' (AZIMUTH  320. 0) ' , 46X, ' ( AZ IMUTH  230.0)',/) 

407  213  FORMAT (/, 2X, ' (AZIMUTH  050. 0) ', 46X, ' (AZ IMUTH  140.0)',/) 

408  300  FORMAT ( 10 (F7. 2, IX) ) 
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7 

409 

400 

F0RMATC9X, \) 

410 

401 

FORMAT (/, 5X, ' 

411 

402 

FORMAT*/, 5X, ' 

412 

403 

FORMAT (/, 5X, ' 

413 

404 

FORMAT (/, 5X, ' 

414 

405 

FORMAT CSX, Al, 

415 

406 

F0RMAT(5X, A5, 

416 

407 

FORMAT (5X, A2, 

417 

406 

FORMAT (5X, A2, 

4ia 

409 

FORMAT (5X, A3, 

419 

410 

FORMAT*/, 5X, ' 

420 

411 

FORMAT (/, 5X, ' 

421 

412 

FORMATC/, 5X, ' 

422 

413 

FORMAT(/, 5X, ' 

423 

414 

FORMATC/, 5X, ' 

424 

415 

FORMAT (/, 5X, ' 

425 

416 

FORMATC/, 5X, ' 

426 

417 

FORMAT (/, 5X, ' 

427 

413 

FORMATC/, 5X, ' 

426 

419 

FORMAT  C/, 5X,  ' 

429 

420 

FORMATC/, 5X, ' 

430 

421 

FORMAT (/, 5X, ' 

431 

451 

FORMAT (/, 5X, ' 

432 

END 
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COLUMNS  OF  OUTPUT  FRAME (N=lto25) :  '  ,\> 

ROWS  OF  OUTPUT  FRAME(M=ito25) :  ',\) 

12, A2, \) 

\) 

12, Al, 12, A3, \) 

13, A3, \) 

\) 

AZIMUTHCO.O  to  360.0  DEGREES):  '  ,\) 

ELEVATION (90. 0  to  -90.0  DEGREES):  ',\) 

NUMBER  OF  SMOOTHINGS:  ',\> 

TRIM (0=N0, l=Xs, 2=Ys) :  '  ,\> 

3,4  OR  3  SUBGRIDS:  '  ,\> 

TITLE  OF  GRAPHCUP  TO  20  CHAR):  ',\) 

DO  YOU  WANT  TO  CHANGE  PARAMETERS'1  '_ ,  \ ) 

DO  YOU  WANT  TO  REPEAT  THE  PROCESS'1  "',  \ ) 

DO  YOU  WANT  FOURIER  TRANSFORMATION  ?  ',\) 

DO  YOU  WANT  TO  MAKE  GRAPH  7  ',\) 

DO  YOU  WANT  CONTOUR  MAP  ?  ',\> 

DO  YOU  WANT  TO  DRAW  CARVE  ?  '  ,  \ ) 

SEND  GRAPH  TO  THE  PRINTERCY  or  N):  ',\) 


Name 


Type 


Offset  P  Class 


ANSWER 

CHAR*1 

AZIM 

REAL 

COS 

CTEXT 

CHAR*20 

DLEV 

REAL 

DM 

REAL 

DN 

REAL 

DY 

REAL 

ELEV 

REAL 

FLOAT 

I 

INTEGER*2 

IDIV 

INTEGER*2 

IMGPAR 

REAL 

I  PRO  J 

INTEGER*2 

ITRIM 

INTEGER*2 

IV 

REAL 

J 

INTEGER*2 

J  J 

INTEGER*2 

K 

INTEGER*2 

L 

INTEGER*2 

LDIG 

INTEGER*2 

LWGT 

INTEGER*2 

M 

INTEGER*2 

MASK 

INTEGER*2 

30 
194 

174 
264 
230 
226 
206 
198 


204 
253 


202 


INTRINSIC 


INTRINSIC 


LARGE 


34 

1  10 

270 

262 

10552 

/WORK 

/ 

10604 

/WORK 

/ 

28 

10656 

/WORK 

/ 
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D  Line**  1 

MM 

INTEGER*2 

N 

INTEGER*2 

NN 

INTEGER*;! 

NRNG 

INTEGER*£ 

p 

REAL 

R 

REAL 

Rl 

REAL 

R2 

REAL 

RLPART 

REAL 

S 

REAL 

sin 

SORT 

TEMP 

REAL 

TRM 

REAL 

U 

REAL 

VERTEX 

REAL 

X 

REAL 

XLOL 

REAL 

XR 

REAL 

XUPR 

REAL 

Y 

REAL 

YLOL 

REAL 

YR 

REAL 

YUPR 

REAL 

Z 

REAL 

ZF 

REAL 

ZFMAX 

REAL 

ZFMIN 

REAL 

ZLEV 

REAL 

ZLQW 

REAL 

M  icros<: 
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£33 

26 

£46 

34 

334 

0 

LARGE 

0 

LARGE 

a 

LARGE 

254 

0 

LARGE 

INTRINSIC 

INTRINSIC 

73 

0 

LARGE 

94 

16656 

/WORK   / 

5403 

/WORK   / 

2 

210 

10 

7323 

/WORK   / 

6 

214 

14 

0 

/WORK   / 

2704 

/WORK   / 

213 

222 

10443 

/WORK   / 

13 

Name 


Type 


Si:e 


Class 


AXIS 

LINE 

MAIN 

MESHS 

NUMBER 

P3D2D 

PLOT 

PLOTS 

SCALE 

STAXIS 

SYMBOL 

WINDOW 

WORK 

ZCNTUR 

ZLEVEL 


16720 


SUBROUTINE 

SUBROUTINE 

PROGRAM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

COMMON 

SUBROUTINE 

SUBROUTINE 
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